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CHAPTER  1 


INTRODUCTION 

The  purpose  of  this  investigation  was  to  extend  present 
understanding  of  aerodynamic  drag  acting  on  cylindrical 
bodies  traveling  concentrically  in  circular  tubes.  The 
primary  motivation  for  this  study  was  that  there  existed  a 
lack  of  detailed  information  related  to  the  near  flow  field 
surrounding  a  vehicle  moving  through  an  incompressible  medium 
within  a  very  long  concentric  tube. 

Considerable  interest  has  been  expressed  within  the  mass 
transit  research  community  (Dept,  of  Transportation)  with 
respect  to  the  near  flow  field  encompassing  a  train  during 
travel  through  a  very  long  tube  (tunnel).  There  is  also 
interest  in  the  U.S.  Navy  concerning  the  mathematically 
similar  problem  of  pushing  a  torpedo  out  of  a  long  over— sized 
tube.  The  work  presented  in  this  study  has  been  tailored  to 
the  range  of  torpedo-tube  combinations  of  most  interest  to 
the  Navy. 

Since  much  of  the  published  literature  on  this  subject  is 
related  to  train-tunnels,  a  brief  discussion  of  that 
application  is  presented.  Several  types  of  high-speed  mass 
transit  ground  transport  systems  have  been  proposed. 
High-speed  vehicles  (usually  trains)  travel  in  enclosed 
guideways  (usually  tunnels)  for  significant  distances.  Since 
the  vehicles  travel  at  sufficiently  low  velocity  it  is 
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reasonable  to  regard  the  air  medium  as  essentially 
incompressible  (Mach  number  less  than  0.2).  The  present 
analysis  is  then  applicable. 

The  physics  of  the  train-tunnel  and  torpedo-tube  are  such 
that  fluid  is  forced  (by  the  naturally  developing  pressure 
drop)  to  flow  down  the  annular  gap  between  the  vehicle  and 
the  guideway  to  fill  the  void  created  at  the  back  of  the 
vehicle.  Motion  of  the  fluid  is  strongly  coupled  to  the 
motion  of  the  vehicle  through  the  conservation  mechanisms  of 
mass  and  momentum.  Specifically,  an  estimate  of  aerodynamic 
drag  on  the  vehicle  is  desired.  Since  total  drag  is  composed 
of  two  separate  parts, i.e.,  shear  drag  and  pressure  drag,  the 
relative  contribution  of  each  is  also  sought.  It  is 
anticipated  that  future  work  will  investigate  methods  of 
reducing  drag  by  understanding  the  predominant  mechanisms 
using  insight  gained  from  this  study. 

Three  main  investigations  were  carried  out  to  address  the 
problem; 

1.  Integral  Formulation,  presented  in  Chap.  2 

2.  Experimental  Formulation,  presented  in  Chap.  3 

3.  Differential  Formulation,  presented  in  Chap.  4. 

Let  us  discuss  these  three  approaches  here  briefly. 
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I NVE STIGATION  #1 


INTEGRAL  FORMULATION 

A  simplified  control  volume  approach  was  initially  taken. 
Drag  results  are  readily  obtainable  but  lack  the  desired 
local  detailed  information  about  the  velocity  and  pressure 
fields  since  this  method  forces  one  to  deal  with  averaged 
quantities.  It  is  certainly  useful,  however,  to  estimate  the 
terminal  velocity  and  tube  length  required  to  conduct  a  drop 
test  using  a  particular  tube-model  pair. 

A  simplified  engineering  model,  using  hydraulic  loss 
coefficients  and  pipe  friction  factors,  was  constructed. 
Since  this  simplified  model  assumed  that  vehicle  drag 
coefficient  was  constant,  i.e.,  not  a  function  of  velocity, 
the  model  was  called  the  constant-C^  model.  The  constant-C^ 
model  was  based  on  Bernoulli  type  pressure  drops  in  the  nose 
region  (very  little  shear  losses),  Kotlow's  [10]  integral 
method  in  the  straight  annular  region,  and  Hoerner's  [1]  loss 
coefficient  model  for  the  tail.  Total  drag  was  obtained  by 
adding  the  independently  calculated  drag  components  of  the 
separate  regions.  The  result  is  a  simplified  engineering 
model  capable  of  reasonably  reproducing  the  drag  predicted  in 
the  more  exact  differential  approach  described  in 
Investigation  #3,  with  an  enormous  savings  of  time  and 
effort . 
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INVESTIGATION  #2 
EXPERIMENTAL  FORMULATION 

Wind  tunnel  testing  is  not  practical  due  to  the 
requirement  that  the  body  move  with  respect  to  the  fixed 
wall.  A  fixed  body  within  the  wind  tunnel  violates  similarity 
of  boundary  conditions.  Experimental  studies  were  conducted 
using  cylindrical  bodies  with  hemispherical  noses  and  conical 
tails  to  verify  the  analytical  work.  A  suitable  drop  test 
apparatus,  using  a  water  filled  tube,  was  designed  and 
fabricated,  similar  to  Hoppe  and  Gouse  [11]  and  later  by 
Nayak  et  al.  [12].  By  closing  the  end  of  the  tube  it  is 
possible  to  simulate  an  infinitely  long  tube.  This  is  because 
in  an  infinitely  long  tube  there  is  no  induced  axial  flow  in 
the  tube  due  to  its  infinite  resistance  to  flow.  A  suitable 
guideway  and  release  mechanism  was  designed  and  constructed. 
Results  were  obtained  for  gap  Reynolds  numbers  (based  on  the 
annular  gap)  up  to  2E4.  High  speed  photography  was  used  to 
construct  a  position  vs.  time  history  for  the  model  from  rest 
to  terminal  velocity.  Terminal  velocity  was  sufficiently  high 
to  ensure  that  turbulent  flow  had  been  established  in  the 
annular  clearance. 
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I N VE S TIGATION  #3 


DIFFERENTIAL  FORMULATION 

The  first  step  was  a  finite  difference  analysis  to  solve 
the  axisymmetric  inviscid  velocity  and  pressure  fields 
surrounding  the  body. 

The  main  analytical  effort  was  to  conduct  a  differential 
analysis  of  the  viscous  flow  field  over  a  range  of  laminar 
and  turbulent  Reynolds  numbers  using  the  finite  difference 
approach.  After  the  velocity  and  pressure  fields  were 
determined  the  data  were  postprocessed  to  obtain  forces  on 
the  surface  of  the  body. 

From  the  known  velocity  field,  local  shear  stresses  on  the 
vehicle  were  obtained  and  integrated  over  the  surface  of  the 
vehicle  to  obtain  the  overall  axial  component  of  the  shear 
forces,  briefly  called  shear  drag.  From  the  pressure  field 
the  axial  component  of  overall  pressure  drag,  also  called 
form  or  profile  drag,  was  obtained.  Total  drag  was  found  by 
adding  the  overall  shear  and  form  drag. 

Several  simplifying  assumptions  were  made  during  this 

study: 

1.  Incompressible  flow 

2.  Quasi-steady  flow  past  the  body 

3.  Infinite  tube  length 

4.  Axisymmetric  flow 

At  terminal  velocity,  the  drag  on  the  body  equals  the  net 
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weight  (weight  in  fluid  minus  weight  of  displaced  water). 
Conversely,  the  differential  analysis,  in  which  the  body  is 
at  the  terminal  velocity,  should  predict  a  drag  force  equal 
to  the  net  weight  from  the  drop-test  experiment  and,  in  fact, 

excellent  agreement  was  obtained. 

Drag  coefficients  for  hydraulically  smooth,  streamlined 
bodies  may  be  extrapolated  with  confidence  for  turbulent 
Reynolds  numbers  of  2E4  up  to  prototype  values  of  5E6  using  a 
power-law  relation. 

A  curve-fit  of  versus  Reynolds  number  has  been  proposed 
in  this  study  for  the  closed  environment  of  a  tube-vehicle 
for  the  range  of  diameter  ratios  and  turbulent  Reynolds 
numbers  of  interest  to  the  Navy. 
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1.1  Background 

The  problem  of  immediate  interest  to  the  U.S.  Navy  is  the 
development  of  an  accurate  engineering  model  capable  of 
predicting  fluid  drag  on  a  torpedo  during  a  "swimout"  launch 
from  a  submarine  torpedo  tube.  Swimout  simply  means  the 
torpedo  engine  is  started  and  the  torpedo  powers  itself  down 
the  tube  and  exits  into  the  ocean  as  shown^  in  figure  la  . 
figure  lb  schematically  shows  a  torpedo  launching  system, 
with  a  pump  providing  an  external  pressure  gradient  to  push 
the  torpedo  out  of  the  submarine.  The  external  pressure 
gradient  problem  was  not  addressed  in  this  study.  This 
investigation  was  concerned  with  the  flow  problem  shown  in 
figure  la  only,  a  simplified  model  of  which  is  presented  in 
figure  2  . 

This  study  features  a  differential  computer  model  that 
incorporates  a  hemispherical  nose  and  a  conical  tail  fixed  to 
the  straight  cylindrical  section  and  thus  complements  and 
considerably  extends  all  previous  work.  In  addition,  this 
study  included  the  design  of  an  experimental  drop  test 
apparatus  that  was  built  and  utilized  to  collect  time  vs. 


1.  Figures  are  collected  at  the  back  of  this  document. 

2.  Later,  in  figure  8  ,  there  is  shown  a  cylinder  of 
diameter,  d,  and  length,  L,  having  a  hemispherical  nose  and 
conical  tail  on  the  centerline  of  a  closed-end  tube  of 
diameter,  D,  and  length,  L^.  Other  relevant  parameters,  such 
as  absolute  vehicle  velocity,  v,  and  average  flow  velocities, 
u  and  w,  are  shown  and  related  in  figure  12  . 
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distance  data  for  actual  scale  models  in  the  geometrical 
ratios  of  interest  to  the  Navy.  Before  progressing  into  the 
details  of  this  analysis,  a  historical  background  and  the 
significance  of  a  number  of  previous  studies  are  presented, 
most  of  which  have  been  on  the  subject  of  trains  traveling 
through  long  tunnels. 

Gouse  [7]  provides  an  excellent  review  of  the  early  German 
work  on  the  train-tunnel  problem.  The  following  discussion  is 
offered  with  credit  to  [7].  One  of  the  earliest  German  papers 
is  by  weismann  [14],  who  assumed  that  for  an  open  tunnel  drag 
would  be  proportional  to  (V^  -  )  where  U„  is  an  induced 

far-  field  velocity  in  the  open  tunnel  (U„-0  for  a  very  long 
tube).  We  now  know,  however,  from  Davidson  [8]  or  Hammitt 
[20]  and  [35],  that  drag  is  proportional  to  (V  -  U„)  .  Prandtl 
and  some  fellow  researchers  also  did  work  on  train-tube 
aerodynamic  drag.  Tollmein  et  al.  published  a  paper  [15]  on 
the  subject  in  1927,  as  did  Danger  [16].  As  discussed  in  [7], 
this  group  placed  model  train  cars  in  a  stationary  wind 
tunnel.  They  made  the  assumption  that  drag  could  be  found: 

1  2 

F ,  -  -  p  a  C ,  (V  -  U,„)  , 

°  2 

where 

F,  »  model  drag  force 
a 

p  *  density  of  fluid  medium 
a  =  model  cross-sectional  area 

C,  »  model  drag  coefficient 
a 

V  ■  model  absolute  velocity 
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U  -  induced  far-field  velocity  in  the  tunnel. 

<D 

They  assumed  was  only  a  function  of  train  and  tunnel 

geometry  and  Reynolds  number  based  on  (V  -  U^) .  They  ignored 
the  fact  that  they  had  violated  the  boundary  condition  of  a 
moving  inner  wall,  which,  as  discussed  in  Sect.  3.2,  is  very 
significant.  They  also  ignored  friction  on  the  walls  of  the 
tunnel  in  the  vicinity  of  the  train,  which  turns  out  to  be  a 
fairly  reasonable  assumption,  but  they  did  not  try  to  estimate 
its  relative  importance.  They  went  on  to  relate  pressure 
drops  using  the  steady  flow  energy  equation  with  pipe  loss 
coefficients  down  the  tunnel,  which  has  subsequently  been 
shown  by  a  number  of  researchers  to  be  acceptable.  Some  of 
these  early  deficiencies  were  addressed  by  Sutter  [4}.  He 
estimated  the  relative  importance  of  wall  friction  in  the 
vicinity  of  the  train  and  found  it  to  be  negligible  for  large 
clearances  and  short  tunnels.  However,  Kotlow  [10]  has  shown 
that  wall  friction  plays  an  important  role  for  small 
clearances  in  both  the  open—  and  closed— end  tunnels.  Sutter 
also  improved  upon  Tollmein's  work  by  including  the  velocity 
in  the  gap,  w,  in  his  analysis.  Sutter  also  did  extensive 
testing  involving  trains  entering  and  exiting  tunnels.  He 
then  generalized  his  simple  loss  coefficient  theory  to 
include  the  transient  as  well  as  the  steady  state  case.  It 
is  interesting  to  note  from  Sutter's  work  that  maximum  drag 
occurs  approximately  when  the  rear  of  the  train  enters  the 
tunnel  and  then  decreases  to  a  steady  value.  For  a  180  meter 
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train  traversing  a  3300  meter  tunnel,  approximately  2200 
meters  of  travel  are  required  to  reach  the  steady  value.  Drag 
was  calculated  by  measurement  of  engine  power  output.  Rolling 
losses  were  predetermined  by  running  the  train  in  open  spaces 
so  that  the  tunnel  effect  could  be  isolated.  Unfortunately, 
Sutter  did  not  run  the  train  into  a  very  long  tunnel  that  was 
completely  blocked  off  on  the  far  end  and  let  it  coast/brake 
to  a  stop  after  reaching  steady  state.  Under  these 
circumstances  steady  state  would  be  reached  very  soon  after 
the  rear  of  the  train  entered  the  tunnel  so  that  good  data 
would  be  easily  and  safely  collected.  There  was  probably  no 
reason  to  perform  such  a  test  since  a  long  tunnel  in  Europe 
is  2  miles  long.  The  present  study  examines  the  limiting 
case  of  zero  induced  flow,  which  remains  to  the  present 
untested  at  full  scale  Re^  of  lo”^.  Scale  model  tests  such  as 
Davidson  [8]  in  1974  have  been  done  up  to  Re^  of  200,000,  and 
in  1972  the  Jet  Propulsion  Laboratory  at  Caltech  in  Pasadena, 
California  [17],  reported  scale  model  results  up  to  Re^=3E6. 

Abramovich  [5]  and  Kara  [6]  also  developed  loss 
coefficient  models  based  on  Sutter's  data.  More  recently 
analytical  work  has  been  done  by  Gouse  et  al.  17],  Davidson 
[8],  Sud  and  Chaddock  [9],  and  Kotlow  [10].  Ref  [7]  and  [8] 
used  basic  integral  techniques  combined  with  loss  coefficient 
models  to  obtain  drag.  The  significance  of  [8]  was  that 
compressibility  effects  were  measured  and  analytically 
modeled.  The  integral  technique  applied  to  the  incompressible 
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case  was  greatly  improved  by  Sud  [9]  through  use  of  an  eddy 
viscosity  model  in  the  sublayer  and  von  Karman's  similarity 
hypothesis  in  the  turbulent  region.  Sud  extended  the  analysis 
to  include  the  developing  region.  Kotlow  [10]  then  extended 
Sud's  analysis  to  include  a  wide  range  of  blockage  ratios  and 
Reynolds  numbers. 

This  investigation  is  essentially  an  extension  of  work 
done  by  Kotlow  [10]  concerning  fluid  dynamic  drag  on  smooth 
cylindrical  bodies  traveling  concentrically  down  very  long 
fluid-filled  tubes.  Figure  3  shows  the  shapes  of  the 
velocity  profile  at  various  points  in  the  flow  field.  Kotlow 
assumed  the  flow  was  steady  and  incompressible  and  that  all 
walls  were  hydraulically  smooth.  This  study  makes  the  same 
assumptions.  Whereas  Kotlow  assumed  that  flow  enters  the 
straight  annular  region  in  a  rectangular  velocity  profile  as 
in  figure  3  ,  this  study  allowed  for  the  simultaneous 
solution  of  the  flow  field  in  all  regions,  producing  a 
realistically  skewed  velocity  profile  at  the  entrance  to  the 
straight  annular  region. 

Experimental  efforts  to  understand  concentric  turbulent 

annular  tube  flow  were  recorded  as  early  as  1907,  by  Becker 

[50].  Quarmby  [18]  was  successful  in  predicting  friction 

factors  for  fully  developed  concentric  turbulent  annular  tube 

flow.  The  theoretical  results  presented  were  for  the  range 

6E3  <  Re-  <  4.5X10^,  for  six  values  of  P,  ranging  from  1.05 

“  p  “ 

to  50,  Re _*V( D-d )/v .  Quarmby  [18]  developed  a  basic  approach 
P 
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that  was  successfully  extended  to  the  case  of  the  inner  wall 
moving  with  respect  to  the  outer  wall.  The  technique  is  to 
use  Deissler's  equation  for  eddy  diffusivity  to  obtain 
velocity  profiles  for  the  sublayer  and  von  Karman's 
similarity  hypothesis  to  obtain  profiles  in  the  turbulent 
layer.  Quarmby  (18)  succeeded  in  obtaining  reasonable 
agreement  with  his  experimental  results.  Wilson  and  Medwell 
[19]  extended  the  Quarmby  [18]  analysis  by  modeling  the 
hydrodynamic  boundary  layer  growth  in  the  entry  region  for 
concentric  turbulent  annular  tube  flow  using  the  momentum 
integral  technique.  The  velocity  profiles  were  determined  by 
use  of  Reichardt's  expression  for  eddy  diffusivity  which  was 
adjusted  for  calculations  near  the  wall  by  implementation  of 
Van  Driest's  [32]  damping  factor.  The  successful  analytical 
models  of  Quarmby  [18]  and  of  Wilson  and  Medwell  [19]  were 
combined  and  extended  by  Sud  and  Chaddock  [9]  to  provide  for 
a  moving  inner  boundary.  The  drag  coefficient  predicted  by 
Sud  and  Chaddock  at  Re^-10®  was  adjusted  using  the  1/6  power 
law  relation  suggested  by  Gouse  [7]  to  obtain  a  value 
corresponding  to  Re^-10^.  The  predicted  coefficient  was  about 
15%  lower  than  the  value  given  by  Davidson  [8],  which  was 
based  on  a  compilation  of  independent  experimental  results 
produced  by  Grittner,  Smith,  JPL,  Hoppe  and  Gouse.  Davidson 
reduced  and  correlated  these  results  to  a  common  plot  of  drag 
coefficient  versus  area  blockage  ratio,  shown  in  figure  4  . 
Sud  and  Chaddock  speculated  that  neglect  of  entrance  effects 
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and  the  use  of  a  simplified  model  for  the  wake  region  in 
their  analysis  may  have  contributed  to  the  discrepancy  between 
their  results  and  the  correlation  of  Davidson, 

In  the  author's  opinion  the  work  of  Kotlow  [10]  is  the 
best  available  analytical  work  done  on  this  subject. 
Therefore,  a  brief  description  of  Kotlow's  analysis  and  the 
computer  algorithm  has  been  included.  Kotlow  [10]  applied  the 
model  of  turbulent  flow  developed  by  Sud  and  Chaddock  [9]  to 
obtain  velocity  profiles  and  pressure  gradients  at  each  axial 
station  down  the  annulus.  There  are  essentially  two  parts  to 
the  Kotlow  analytical  model,  a  fully  developed  region  and  an 
entry  region.  The  fully  developed  model  will  be  discussed 
first. 

In  the  sublayer,  Deissler ’ s equation  for  eddy  diffusivity 
was  used  to  produce  a  first  order  ordinary  differential 
ecjuation  in  du'^'/dy'*' ,  including  a  damping  factor,  n  ,  obtained 
from  concentric  tube  experiments  that  accounts  for  the  effect 
of  the  wall  on  turbulence.  Runge-Kutta  numerical  integration 
was  used  to  obtain  the  dimensionless  velocity  profiles  in  the 
sublayer.  In  the  turbulent  layer  model  Kotlow  (10]  used  a 
second  order  ordinary  differential  equation  in  d^u''’/dy''’^ , 
obtained  by  Von  Karman's  similarity  hypothesis  combined  with 
the  Prandtl  mixing  length.  Runge-Kutta  was  employed  to 
perform  the  integration  and  produce  dimensionless  velocity 
profiles.  Integration  is  continued  as  far  as  the  axis  of  zero 
shear.  Once  the  dimensionless  velocity  profile  was  found  in 
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the  fully  developed  region,  local  vehicle  shear  stress  is 
calculated.  Profile  or  pressure  drag  was  found  using  a  force 
balance  on  the  fluid  in  the  annular  gap. 

The  momentum  integral  technique  was  used  to  locate  each 
profile  in  the  entry  region,  producing  shear  stress  as  a 
function  of  axial  position.  Entry  shear  drag  was  found  by 
integrating  local  stresses  over  the  surface  area.  Profile 
drag  in  the  entry  region  was  found  by  determining  the  local 
pressure  gradient,  computed  by  assuming  irrotational  flow  in 
the  inviscid  core  region,  and  assuming  zero  radial  pressure 
gradient. 

sirs  several  important  advantages  of  the  differential 
approach  used  in  this  study  over  the  integral  techniques 
employed  by  Kotlow  [10].  The  main  advantage  is  that  by 
solving  the  entire  flow  field  around  the  body  in  a  single 
operation  it  is  not  necessary  to  make  many  of  the  assumptions 
necessary  to  carry  out  an  integral  analysis.  For  example,  the 
velocity  profile  at  the  entrance  to  the  straight  annular 
region  was  assumed  by  Kotlow  to  be  "slug"  flow,  figure  (3). 
This  study  showed  that  the  velocity  profile  at  this  location 
is  actually  skewed  due  to  shear  effects  from  the  nose.  This 
skewing  produces  the  effect  of  shortening  the  developing  or 
entry  length.  The  differential  calculations  done  in  this 
study  also  do  not  require  the  assumption  of  irrotational  flow 
in  the  inviscid  core  in  the  entry  region,  in  order  to 
calculate  the  pressure  gradient.  There  are  many  other 
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assumptions  that  can  be  avoided  in  the  matching  and  checking 
of  velocity  profiles,  the  details  of  which  are  very  involved. 

In  the  wake  region,  a  pressure  coefficient  was  calculated 
by  Kotlow  [10]  similar  to  Sud  and  Chaddock  [9]  and  Hammitt 
[20]  based  on  Hoerner's  data  [1]  of  a  model  in  free  field 
flow.  Wake  drag  was  calculated  by  multiplying  this  pressure 
drop  (rise)  by  the  projected  vehicle  area.  This  study  uses  a 
sudden— expansion— in— pipe-flow  analogy  of  the  wake,  in  the 
integral  analysis  in  Chap.  2.  The  differential  solution  in 
Chap.  4  of  this  study  has  the  distinct  advantage  of  not 
requiring  the  use  of  open  environment  data,  or  the  pipe— flow 
analogy. 

Kotlow' s  entry  analysis  showed  that  wall  shear  stresses 
start  high  at  the  inlet  and  rapidly  decay  in  the  axial 
direction.  The  developing  pressure  gradient  also  is  quite 
high  at  the  inlet  but  decays  more  slowly  than  shear  stress. 
Consequently  profile  drag  is  more  significant  than  shear  drag 
in  the  entry  length.  Kotlow  concludes  that,  quantitatively, 
the  entry  length  is  nearly  negligible,  since  it  raises  the 
total  drag  by  an  inconsequential  3.0%.  Based  on  this 
important  conclusion  it  was  assumed,  in  Chap.  2,  that  the 
flow  was  fully  developed  in  the  entire  annular  region. 

Between  1920  and  1960  there  was  essentially  no  progress 
made  until  interest  in  smaller  clearances  and  higher  speeds 
started  to  arise.  Much  experimental  work  has  been  done, 
especially  by  the  Japanese,  who  are  interested  in  high  speed 
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streamlined  trains.  Hoppe  and  Gouse  [11]  provide  an  excellent 
summary  of  early  Japanese,  French  and  Russian  work.  The 
following  discussion  is  provided  with  credit  to  [11].  In  the 
first  important  Japanese  analysis,  Kawaguti  [21]  used  the 
Rayleigh  analogy  in  an  incompressible  laminar  flow  analysis 
to  obtain  the  drag  of  a  train  concentric  to  the  tunnel. 
Kawaguti  then  extended  the  laminar  analysis  to  the  turbulent 
case  by  the  use  of  the  eddy  viscosity  concept.  This  general 
approach  was  later  followed  by  Gouse  [2].  Kawaguti  used  the 
experimental  data  of  Miki  [22]  to  obtain  the  proper  eddy 
diffusivity  coefficient.  Kara  [23]  analytically  predicted  the 
aerodynamic  forces  acting  on  a  train  entering  a  tunnel. 
Actual  drag  data  presented  by  Kawamura  and  Ono  [24]  had  been 
previously  taken.  A  transient  one-dimensional  inviscid  flow 
problem  taking  into  account  compressibility  was  solved.  A 
pressure  wave  propagates  to  the  opposite  end  of  the  tunnel 
(and  is  then  reflected  back)  when  the  nose  of  the  train 
enters  the  tunnel.  The  solution  by  Kara  continues  until 
either  the  reflected  pressure  wave  comes  back  to  the  nose  of 
the  train,  or  the  rear  of  the  train  has  entered  the  tunnel. 
By  assuming  a  value  for  the  train  skin  friction  coefficient, 
taken  to  be  constant  and  equal  to  the  free  environment  value, 
the  pressure  drag  coefficient  and  pressure  distribution  along 
the  train  axis  was  determined.  Pressure  data  were  taken  by 
Kara  and  Okushi  [25]  while  pulling  a  cylinder  (d=  0.63  in, 
L=  2  ft)  through  a  tube  (D=  1.1  in,  L^=  19  ft).  The  model  was 
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suspended  on  a  piano  wire  that  was  stretched  down  the 
centerline  of  the  tube.  Measurements  of  static  pressure  along 
the  wall  of  the  tube  were  taken  and  correlated  with  various 
train  nose  and  tunnel  entrance  shapes.  Having  the  pressures 
along  the  wall  allows  for  direct  prediction  of  profile  drag. 
Hara's  original  analytical  model  underpredicted  profile  drag 
since  it  was  based  on  the  lower  value  of  free  environment 
skin  friction  coefficient. 

Kara  [6]  set  out  to  improve  the  original  analytical  model 
by  getting  a  better  estimate  of  the  transient  skin  friction 
coefficient  as  the  model  enters  the  tube.  To  accomplish  this 
the  pressure  data  mentioned  above  were  used  together  with  a 
one— dimensional  analysis  of  the  annular  flow  between  the 
model  and  the  tube.  This  method  turned  out  to  be  very 
successful  and  Kara  [26]  presented  pressure  data  and  skin 
friction  coefficients  for  specific  Japanese  railroad 
train-tunnel  combinations.  To  complete  his  drag  model  Hara 
[27]  incorporated  the  results  of  Miki  [22]  and  Tanaka  [28] 
for  wind  tunnel  measurements  of  train  afterbody  drag.  Hara 
[27]  then  offers  an  empirical  drag  coefficient  equation, 
developed  for  the  unique  Japanese  geometry  (0=1.75). 

Another  Japanese  author  of  note  is  Fukuchi  [29] ,  who 
studied  drag  on  trains  in  very  long  tunnels.  The 
one— dimensional  approach  of  Hara  was  utilized.  A  drag 
coefficient  model  of  reasonable  accuracy  was  presented.  The 
difficulty  was  in  getting  data  for  very  long  tubes,  so  the 


model,  being  inherently  empirical,  was  less  accurate  than 
Hara's  model.  Fukuchi  [30]  used  the  Hara  approach  to 
calculate  drag  on  trains  passing  through  underground  train 
stations,  and  for  the  case  of  trains  traveling  in  opposite 
directions  inside  the  tunnel. 

Similar  one— dimensional  models  have  been  made  for  the 
French  railroad  by  Nouvion  [31],  and  for  the  Russian  railroad 
by  Abramovich  [32].  All  reported  results  are  empirical  and 
highly  specific  to  the  actual  train— tunnel  being  tested. 

Cromack  [33]  obtained  data  on  1/24^^  scale  slot  cars  in 
the  open  and  going  through  tunnels  in  air.  He  was  able  to 
measure  the  instantaneous  electrical  power,  P,  required  by 
the  car  and  he  obtained  instantaneous  velocity,  V. 
Aerodynamic  drag,  F^,  was  then  calculated  from,  F^  =  P/V.  The 
result  is  .  good  qualitative  trends  between  different 
configurations  but  quite  low  quantitative  predictions. 
Perhaps  the  fault  is  that  the  physical  relation  of  the  car  to 
the  track  is  too  eccentric  to  be  expected  to  conform  to  the 
standard  concentric  predictions.  Eccentric  tubes  tend  to 
require  less  pressure  drop  to  produce  the  same  flow  rate 
through  the  annulus. 

Experiments  done  by  Gregorek  and  Engle  [34]  involve  the 
deceleration  of  cylindrical  models  in  a  long  horizontal  tube. 
The  models  were  propelled  into  the  tube  and  the  deceleration 
was  measured.  Drag  force  is  then  inferred  from  the 
deceleration  of  the  model.  Three  model  diameters  were  tested 
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(d-1.806  in,  d-  2.089  in,  d=  2.414  in)  and  three  L/d  ratios 
(10,20,30).  The  tube  (D=  2.8  in,  L^=  180  ft)  produced  8-1.55, 
1.34,  1.16  respectively.  The  models  all  had  a  streamlined 
boat  tail  afterbody  and  either  a  flat  or  ellipsoidal  nose. 
The  results  of  this  experiment  agree  with  the  generally 
accepted  results  shown  in  figure  4  . 

Hammitt  [35]  presents  a  basic  one-dimensional  analysis  for 
the  induced  flow  produced  by  a  train  traveling  through  a 
tunnel  for  both  the  steady  and  accelerating  cases.  He  assumes 
that  the  vehicle  drag  coefficient  is  constant,  i.e., 
unaffected  by  changing  vehicle  velocity.  This  assumption 
neglects  the  important  Reynold's  number  dependence,  which  most 
investigators  take  as  an  axiom. 

Perhaps  the  most  impressive  experimental  effort  to  date 
as  discussed  by  Kurtz  [17]  was  the  VICS  (Vehicles  in  Confined 
Spaces)  study  conducted  by  Caltech,  JPL  in  1971—1973,  for  the 
DOT  (Dept  of  Transportation).  The  ViCS  facility  consists  of 
three  elaborately  instrumented  long  vertical  tubes,  all  D= 
2.0  in,  the  VICS-32,  VICS-70  and  VICS-120,  where  the  number 
denotes  tube  height,  L,  in  feet.  The  VICS— 70  is  capable  of 
producing  Re^  up  to  10^.  The  investigators  at  JPL  recognized 
that  the  actual  near  flow  field  is  far  more  complicated  than 
the  usually  employed  one-dimensional  analysis  allows.  Their 
idea  was  to  develop  a  long  vertical  test  facility  so  that 
terminal  velocity  could  be  measured  with  a  high  degree  of 
precision.  They  also  wanted  detailed  pressure  data.  In  the 
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VICS-70  models  are  launched  from  a  spring  powered  short 
stroke  catapult  accelerator  to  a  velocity  close  to  the 
terminal  velocity.  The  first  50  feet  of  travel  allows  for 
small  variations  in  model  velocity  to  disappear  and  the  model 
to  settle  down  to  the  terminal  velocity.  The  final  20  feet  of 
travel  are  where  most  of  the  wall  pressure  and  velocity 
readings  are  taken,  in  order  to  maintain  pressure  integrity 
in  the  tube  during  pressurized  runs,  a  system  of  copper  wire 
coils  was  placed  around  the  outside  of  the  tube.  Permanent 
magnets  were  placed  in  the  model  in  order  to  induce  a  voltage 
in  the  coil  as  the  model  passed  by.  Velocity  was  found  from 
V*=  (distance  between  coils)/time.  Static  wall  pressure  was 
found  from  very  accurate  pressure  transducers  capable  of 
giving  a  pressure  profile  as  the  model  passed  by.  Test 
Reynolds  numbers  were  increased  by  using  CO2  and  Freon-12. 
More  information  on  this  method  of  increasing  Re^  is 
contained  in  Appendix  A.  The  facility  was  hydrostatically 
tested  to  300  psig  but  operated  only  up  to  165  psia  (11  atm). 
Use  of  Freon  at  70  psig  and  70®F  is  delicate  because  it  will 
liquefy,  so  a  toggle  valve  was  installed  at  the  bottom  of  the 
tube  to  check  for  this  possibility.  Models  typically  weighed 
25  pounds  and  were  released  with  no  spin. Their  fall  was  guided 
by  molded  plastic  runners  to  remain  concentric  (or  eccentric, 
if  desired).  A  special  arrestor  was  constructed  of  cellulose 
sponge  sections  to  absorb  the  model's  kinetic  energy  at  the 
end  of  the  run. 
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In  order  to  achieve  Re^  =  3E6,  for  incompressible  flow, 
the  VICS-120  was  constructed.  It  is  attached  to  the  outside 
of  the  9  story  Space  Sciences  Building  (Bldg.  183)  at  the 
JPL.  It  is  in  a  class  all  by  itself,  even  considering  the 
impressive  VICS-70  just  described.  It  is  capable  of  operating 
at  50  atm  and  uses  CO2,  N2,  Freon-12  or  air  as  the  working 
fluid.  It  is  very  flexible  in  terms  of  being  able  to  test 
configurations  requiring  simultaneous  movement  of  models  in 
interconnecting  tubes,  as  would  be  experienced  in  underground 
train  systems.  Porous  wall  sections  can  be  added,  as  can  wall 

elements  with  different  roughness. 

A  gas  gun  catapult  can  launch  a  50  pound  model  at  100  fps 
in  3  ft,  with  no  spin.  The  model '  s  travel  was  guided  by  plastic 
runners  down  the  centerline  of  the  tube  (D*  2",  120 
ft).  Models  typically  travel  at  10  -  100  fps.  Velocity  is 
measured  by  optical  sensors  used  to  form  what  Kurtz  [17] 
calls  speed  traps  created  by  the  model  breaking  a  light 
beam.  The  speed  traps  have  a  time  resolution  of  0.10  ms  and 
an  accuracy  of  0.16  ms  which  allows  precise  velocity 
measurement  up  to  250  fps.  Pressure  signatures  are  produced 
by  fast  response  (200  Hz)  pressure  transducers  mounted  on  the 
wall  of  the  tube.  Great  ingenuity  went  into  the  design  of 
models  intended  to  simulate  actual  train  cars.  Distributed 
roughness  produced  by  coating  the  exterior  of  the  model  in 
0.003  diameter  polystyrene  balls  was  used.  Special  grooving 
was  fashioned  to  simulate  door,  windows  and  other  openings. 
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The  results  of  the  JPL  tests  for  the  conditions  assumed  in 
the  present  study  of  infinitely  long  (  closed-end  ) tubes  are 
compared  to  the  results  of  the  differential  section.  Chap.  4, 
of  this  study.  Good  agreement  was  realized. 

The  author  considers  that  many  of  the  design  concepts  that 
went  into  the  VICS-70  and  VICS-120  are  very  significant  with 
respect  to  the  topic  under  study.  The  author  researched  most 
of  the  parameters  that  control  the  effectiveness  of  any 
experimental  effort  related  to  this  topic.  Significant  topics 
include:  choice  of  fluid  medium,  model  design,  guideways, 
release  mechanisms,  arrestors,  velocity  measurement  and  model 
stability  considerations.  The  author  of  this  report  reached 
the  conclusion  that  the  JPL  design  is  optimum  in  all  these 
respects  for  steady  state  testing;  however,  it  does  have  some 
important  limitations  which  are  briefly  discussed. 

The  main  objective  of  this  study  was  to  perform  a  detailed 
study  of  the  flow  field  surrounding  the  model;  this  cannot  be 
done  by  JPL,  or  any  other  test  facility.  The  secondary 
objective  of  this  study  was  to  drop  test  models  having 
specific  geometries  and  values  of  dimensionless  parameters  to 
determine  terminal  velocity  and  drag.  The  drop  test  apparatus 
used  in  this  study.  Section  3.3,  produced  accurate, 
repeatable  terminal  velocity  data,  as  verified  by  the 
differential  simulation.  Section  4.8.  Since  the  Navy  is 
interested  in  the  launching  of  torpedoes,  it  is  natural  to  be 
interested  in  transient  velocity  data.  In  this  respect,  the 
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drop  test  done  in  this  study  had  the  advantage  over  the 
steady-state-only-JPL  in  that  transient  model  velocity  data 
werecollected  and  are  presented  in  Section  3.10. 

The  JPL  facilities,  while  capable  of  providing  good  tube 
wall  pressure  data,  are  not  capable  of  providing  pressure 
data  at  the  model  surface.  In  addition,  JPL  is  incapable  of 
directly  providing  detailed  surface  stress  data  on  the  model. 
The  best  they  can  do  is  back-calculate  a  total  shear  drag  by 
subtracting  pressure  drag  (assuming  zero  radial  pressure 
gradient  and  using  wall  pressures)  from  total  drag.  Another 
major  area  that  cannot  be  confidently  addressed  by  JPL  is  the 
question  of  separation  at  the  rear  of  the  model. 

The  differential  simulation  done  in  this  study  overcomes 
all  of  the  limitations  just  discussed,  except  the  transient 
case,  which  was  addressed  experimentally.  The  author  suggests 
that,  by  providing  information  on  local  surface  distributions 
of  pressure  and  shear  as  well  as  knowledge  of  separation, 
this  study  significantly  contributes  to  the  overall  body  of 
knowledge . 
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CHAPTER  2 

CONTROL  VOLUME  FORMULATION 
Sect.  2.1  Basic  Assumptions 

The  purpose  of  this  chapter  is  to  present  the  first  of 
three  studies,  a  control  volume  formulation  of  the  basic 
equation  of  motion  for  a  cylindrical  body  falling  through 
fluid  in  a  long  tube.  This  differential  equation  will  then  be 
solved  both  analytically  and  numerically  and  subsequently 
compared  with  both  the  experimental  and  numerical 
(differential)  formulations. 

The  basic  assumptions  mentioned  in  Chap.  1  will  be  adopted 
here : 

1.  The  flow  is  incompressible  and  viscous. 

2.  All  surfaces  are  hydraulically  smooth. 

3.  Unsteady  hydrodynamic  effects  are  modeled  by  a 
hydrodynamic  mass  term,  m^^,  developed  in  Appendix  D. 

4.  The  body  falls  vertically  and  concentrically. 

The  shear  and  pressure  forces  on  the  body  will  be 
approximated  by  simple  one— dimensional  flow  arguments. 

Sect.  2.2  Application  of  Newton's  Law 

A  schematic  of  the  geometry  of  the  falling  body  is  shown 
in  figure  12  .  Newton's  second  law,  figure  14  ,  states  that 
the  summation  of  vertical  forces  on  the  body  is  proportional 
to  its  vertical  acceleration: 
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2-1 


Z  F  - 


(m  + 


where  W„  -  net  weight  -  W  - 

F.  -  drag  force  -  see  eq.  2-2 
d 

F  -  buoyant  force  =  pgA 
b 


Id  *  body  mass  =  W/g 

=  hydrodynamic  mass 
y  B  body  velocity. 


=  see  eq.  2-4c 


The  body  drag  force  can  be  modeled  in  terms  of  an  overall 
coefficient  referred  to  body  velocity  and 

cross-sectional  area: 


P  2 

^d  -  ^d  “  - ■ 

in  general,  for  a  smooth  surface,  is  a  functi 

shape,  Eeynolds  number,  and  tube-to-body  diameter 

where  S-D/d.  differential  equation 

Combining  eq.  2-1  ana 

for  body  motion  becomes 


dV 

.  i  p  a  Cd  v"  t  (m  *  mj,)  - 


.dh.ect  to  the  initial  condition  V-O  at  t-O.  Fq.  2-3  may  be 

solved  analytically,  if  C,  is  assumed  constant,  a 
numerically,  if  C,  varies  with  Eeynolds  —  - 
aolutions  may  then  be  used  to  fix  design  conditions 

experimental  formulation. 
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Sect.  2.3  Solution  for  Constant  Drag  Coefficient 

Assume  that  the  drag  coefficient,  C^,  is  constant  and 
equal  to  the  drag  coefficient,  C^^.,  at  terminal  velocity,  V^. 
Eq.  2-3  may  then  be  solved  in  closed  form  by  separation  of 
variables  and  integration: 


V  «  tanh(bQt). 

The  parameter  is  the  terminal  velocity  of  the  falling 

body,  which  is  reached  within  99.5%  when  >  3. 


Vi.  “  t>„ 

t  o  o 


2W. 


n 


LCdt  P  aj 


1/2 


and 


2  ( m  +  mj^ ) 

P  ^dt  ^ 


(Wn  P 


1/2 


2-4a 


2-4b 


m  +  m. 


a  value  for  m^^  may  be  found  from  Appendix  D: 

P  A 


m. 


2  ' 
r-1 


where  A  -  volume  of  the  model . 

The  displacement  and  acceleration  are  obtained  from 
S  -  JVdt  =  a^  In  [cosh(b^t)]. 


dV 

and  —  =  a 

dt 


2-4c 


eq.  2-4 


2-5 


Lcosh(bQt ) J 
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wi  rioQian  estimates  of  body 
2  4  and  2-5  provide  reasonable  desig 
E,.  2-4  and  P  action,  which  develops  an 

.otion  if  combined  with  the  next  . _„„body 


.roximate  method 


for  predicting  the  (nearly  constant)  body 


drag  coe 


fficient,  Cj. 


sect  0  4  nrao  Coefficient  Estima^  ,,45 

oross-sectrons  2-3-4-5, 

body  flow  field  into  three  regions:  the 

which  divide  the  near-  y  and  the  tall 

nose  region  (2-3).  the  straight  sectron  (3  ). 

region  (4-5).  coefficient  as 

1-he  overall  drag 

It  is  proposed  to  model  the  ove 

.  1-ail  contributions: 

raf  nose  straight-section,  and  tai 

the  sum  of  nose,  bULo  ^ 


=n  *  "^s  '^t 


... ...... .........  “• "" 

one-dimensional  and  steady,  in  a  frame  relative 

ib-s  frame  the  fluid  approaches  section  2 

Tellcitr  V  and  flows  into  the  straight  (gap)  section  3-4,  at 

velocity,  W,  called  gap  velocity,  given  by  continuity: 


W  -  V/Yf 


where 


A-a 

.  tube  area 


28 


Sect.  2.5  Nose  Drag  Contribution 

Write  the  one-dimensional  steady  flow  energy  equation 
across  the  nose  section  (2-3): 


where  HL  denotes  the  head  loss.  We  approximate  section  2-3 
as  a  'duct  entrance'  and  use  the  well-rounded  entrance  head 
loss  f rom,  e. g. , White  [39]: 


HL 


2-3 


w 

K„  —  (1  -  y)^  »  K  =  0.05 

n  2  n 


2-8 


We  neglect  shear  stress  on  the  body  nose  and  assume  that  the 
nose  drag  is  given  by  =  (p^  -  where 

1  +  (1-Y)^  ■ 

2 

Y  J 


This  is  equivalent  to  a  nose  drag 


coefficient  of 


C 


n 


2(Pn-P3) 


1 

~2 

Y 


+  ^n 


(1-y)^ 


2 

Y 


2-9 


For  example,  when  P^D/d^l.S,  eq.  2-7  gives  y“0.6094,  and  then 
eq.  2-9  gives  Cj^=2.71. 
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Sect.  2.6  Tail  Drag 

Write  the  steady  flow  energy  equation  across  the  tail 
section  ( 4-5 ) : 


P 


+ 


£5 

P 


+ 


+  HL4_5 


For  the  tail  section  we  may  approximate  the  head  loss  as  a 
'sudden  expansion'  where,  from  White  [39]: 


HL4_5 


K^.  =  1.0  . 


2-10 


In  a  similar  manner  to  the  nose-drag  estimate,  we  approximate 
the  drag  on  the  tail  as  =  (p^  -  P5)a»  where 


(P4  -  P5) 


(y^-1)  +  (1-y)^ 


2 

Y 


This  is  equivalent  to  a  tail  drag  coefficient  of 


2  (P4-P5) 
P 


2  1 
Y  -1 


(1-y) 


+  K, 


2-11 


Continuing  our  example  of  3=1.6,  y  =  0.6094,  then  eq.  2-14 
predicts  that  =  -1.28,  where  the  negative  sign  indicates  a 
thrust  or  forward  force. 


Sect.  2.7  Straight  Section  Drag  Contribution 


A  control  volume  enclosing  the  fluid  within  the  gap 


V 


r 


section  (3-4)  in  figure  12  shows  that  there  is  no  average 
momentum  change  because  the  entrance  and  exit  velocities  are 
equal,  v^-v^-w.  Then  the  net  pressure  force  on  the  gap  fluid 
is  balanced  by  shear  forces  on  the  inner  and  outer  surfaces 
of  the  gap  fluid; 


z>,  z 

(P3  -  P4)(A  -  a)  -  ltd 

This  pressure  drop,  induced  by  flow  through  the  gap,  causes  a 
pressure  drag  on  the  body,  Fp  *  (p^  —  P4)a.  In  addition, 

there  is  a  shear  force  on  the  body,  ■«  itd/x^dz,  caused  by 

the  inner  shear  stress,  .  The  total  drag  on  the  straight 
section  is  then  f  *  F  +  F  . 

5  p  c 


r  ^ 


n 


dz  + 


llD 


dz  . 
o 


2-12 


."4 


=  (Po  -  P/i)a  + 


dz  . 


2-13 


To  complete  the  analysis  of  the  straight  section  requires 
correlations  for  (P3-P4)  and  .  A  simple  approximation  would 
be  to  form  a  Reynolds  number  based  on  gap  velocity,  w,  and 
use  the  Moody  chart,  e.g..  White  [39],  to  estimate  both 
pressure  drop  and  wall  shear.  However,  Kotlow  [10]  has  made 
an  extensive  study  of  this  gap  flow,  using  a  developing— flow 
analysis  combined  with  a  turbulent  eddy  viscosity 
formulation.  His  results  are  thought  to  be  the  most  accurate 
available  in  the  literature. 
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Sect.  2.8  The  Gap-Flow  Analysis  of  Kotlow 

First,  Kotlow  [10]  determines  that  little  error  is 
incurred  by  taking  equal  to  a  constant.  Thus,  the  integral 
in  eq.  2-13  is  replaced  by  (ndXj^Lg),  where  is  the  length 
of  the  straight  section,  figure  12  .  He  then  gives 
correlations  involving  two  different  Reynolds  numbers: 


1.  For  the  body: 

2.  For  the  gap  flow:  Re 


Re^j  =  Vd/v, 


wD 


yv  -  p^Re^/(p+l) ,  2-14 


where  «  (D-d)  is  the  hydraulic  diameter  of  the  gap.  Kotlow 


correlates  the  pressure  drop  in  the  style  of  a  Moody  friction 
factor : 


P3  -  P4 


’"s  2 

—  -  f^  , 

D  2 
e 


2-15 


where 


3^1 


L2(3^-1) 


2  . 


2.8 


0.001  + 


3.1 


(logioReg) 

Meanwhile,  Kotlow  correlates  the  inner  wall  shear  stress  in 
terms  of  both  body  and  gap  Reynolds  numbers,  as  follows: 


C(3+l) 


T  .  = 

1 


-  f  _  v‘ 

4(3+^;)  ®  2 


2-16 


where 


1.005  3^  , 


1.477  +  0.2195  log^pRe^ 


7,2 


and 


Kotlow  [10]  states  that  these  correlations  are  accurate  to 
+3%  in  the  range  10^<Re^<10®  and  1.0<p<2.0.  Once  (P3-P4)  and 
T.  are  computed  from  eq.  2-15  and  2-16,  the  straight  section 
drag  coefficient  is,  by  definition. 


2  F, 


2-17 


=>  ^  IT* 

a  p  V 

These  correlations  are  somewhat  complex,  but  their  evaluation 
is  straightforward. 

Sect.  2.9  A  Numerical  Example 

Continuing  our  example  of  S“1.6,  or  d/D*0.625,  let  us 
assume,  guided  by  experimental  considerations,  a  body  of 
diameter  d*0.625  in.  moving  in  a  tube  of  diameter  D-1.0  in. 
at  a  velocity  of  56  in/s  or  4.67  ft/s,  in  water  with 
p=1.94slug/ft^  and  v«l . 07E-05f t^/s .  Then  the  gap  velocity,  w, 
is  7.67ft/s,  whence  Re^-22700  and,  since  0^-0.375  in.. 
Re  *=22400.  Let  the  straight  section  length  be  L  =5.43  in. 

’  2 

From  eq.  2—15,  fg=0.0396,  whence  (p^- P4)=12.1  Ibf/ft  . 

2 

From  eq.  2-15,  J;-3.15,  whence  t^=0.361  Ib^/ft  .  Then  the  drag 
on  the  straight  section  is,  from  eq.  2-13, 


n  0.625^ 

(12.1) - + 

4  144 


0.625  5.43 

n  -  (0.361)  - 


12  12 

-  0.0258  +  0.0267  «  0.053  Ib^  .  2-18 

This  is  equivalent,  from  eq.  2-17,  to  a  gap-section  drag 

coefficient  of  C  *1.17.  By  comparison,  a  calculation  using 

s 
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these  same  input  values  in  Kotlow's  eddy-viscosity  computer 
code  gives  Fg=0,058  Ib^  and  €^=1.29/  or  10%  higher.  The 
primary  reason  for  the  error  is  that  the  body  Reynolds  number 
of  22700  is  below  the  value  of  10^  specified  by  Kotlow,  whose 
correlations  were  intended  for  large-scale  bodies. 


Sect.  2.10  Overall  Drag  Coefficient 

Combining  the  nose,  tail,  and  gap  sections,  we  estimate 
the  overall  drag  coefficient  of  the  body  from  eq.  2-9  and 
2-11  to  be 


^n  ^  ‘^s  "  ^t  ' 

+  (Kj^+K,.  )  (1-y) 


+  C. 


2-19 


For  the  models  tested  here,  arid  K^=l,  so: 

Y^  +  (1-y)^ 

”  2  ' 

Y 

where  must  be  computed  from  Kotlow' s  correlations 
to  2-17. 


2-20 


eq.  2-15 


For  our  continuing  example  with  P=1.6,  eq.  2-20  predicts 
that  the  overall  C^=2.58.  We  will  compare  this  estimate  in 
Chap.  3  with  the  experimental  results. 


Sect.  2.11  Application  to  Design  of  the  Experiment 

The  running  numerical  example  used  here  approximates  a 
body-tube  combination  proposed  for  experimental  formulation: 
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d-0.625  in.,  D-1 . 0  in.,  L  -5.43  in.,  made  of  aluminum,  with 
W-0.164  Ibj  and  volume  A-1.64  in^.  We  wish  this  body,  when 
dropped  in  water,  to  reach  98%  of  terminal  velocity, 
V^=56in/s,  in  a  distance  of  about  2  ft  -  to  stay  within  the 
field  of  view  of  the  high  speed  movie  camera. 

We  have  estimated  above  for  this  body  that  C^=2.58  at 
V-56in/s.  We  now  use  this  value  in  the 

constant-drag-coefficient  formulas,  eq.  2-4  and  2-5,  to  check 
our  estimates.  The  net  weight  of  the  body  in  water  is: 

W  -  W  -  pgA  -  0.164  -  (1.94){32.2)(1. 64/1728)  -  0.105  lb,, 
n 

The  hydrodynamic  mass,  m^^,  is  estimated  from  eq.  2-4c,  the 
formula  developed  in  Appendix  D: 

mj^  =  pA/(8^-l)  =  (1.94)  (1.64/1728)/(1. 6^-1)  -  0.00118  slug, 
or 

m  +  mj^  -  0.164/32.2  +  0.00118  =  0.00627  slug. 

2 

The  body  cross-section  area  is  0.00213  ft  .  Eq.  2-4  may  be 
used  to  predict  body  terminal  velocity: 

V^-  {2(0. 105)/[(2.58)(1.94)(0. 00213)]}^-^^  =  4.44ft/s  -53in/s. 

This  estimate  is  about  5%  less  than  our  initially  guessed  value 
of  56  in/s.  If  necessary,  one  could  iterate  to  find  a  new 
Reynolds  number  and  repeat  the  procedure.  Or  one  could  modify 
body  weight, etc.  ,  to  meet  these  conditions. 
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From  eq.  2-4,  the  parameters  a^=1.17  ft  and  b^=3.78  s 
the  proposed-design  input  values.  We  see  from  2-4  that  V 
-  0.98  V^.  when  (b^  t)  -  2.298;  that  i  s,  tanh(  2 . 298 )  -  0.98. 

Then,  from  eq.  2-5,  the  distance  traveled  to  reach  this 

velocity  is 

S  *  a  ln[cosh(b  t)]  =  1.17  ln[  cosh(  2 . 298 )  ]  =  1.89  ft  *=  2_3 _ iji. 

o  o 

This  is  within  the  2— ft  field  of  view  of  the  camera. 
Therefore,  the  proposed  design  seems  adequate,  at  least  by 
this  simple  constant-C^^  analysis. 

This  estimate  of  distance  traveled  should  actually  be 
conservatively  large,  since  we  have  used  at 

terminal  velocity,  whereas  the  actual  C^(t)  begins  higher 
than  and  continually  decreases  to 

Sect.  2.12  Solution  for  Variable  Drag  Coefficient 

Since  Kotlow's  gap-section  correlation  predicts  that 
varies  with  Reynolds  number,  it  is  expected  that  more 
accuracy  could  be  obtained  by  solving  eq.  2—3  numerically 
using  a  variable  as  predicted  by  eq.  2-15  to  2-20.  A 

computer  program  that  uses  the  Runge— Kutta  method  to 
integrate  eq.  2-3  is  given  in  Appendix  E.l.  It  is  readily 
adaptable  for  any  tube/model  configuration  that  is  desired. 
This  program  was  applied  to  the  proposed  design  body— tube 
combination  discussed  above  and  will  be  compared  in  Chap.  3 
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with  the  experimental  data  and  with  the  constant 
approximation . 

Summary 

This  chapter  has  presented  a  control  volume  analysis  of  a 
cylindrical  body  falling  concentrically  in  a  vertical  tube. 
Using  one-dimensional  flow  approximations  plus  a 
turbulent-gap-flow  study  by  Kotlow  [10],  simple  formulas  are 
proposed  for  overall  body  drag,  and  a  simple  solution  is 
given  for  velocity  and  displacement  of  the  body.  These 
formulas  are  used  to  verify  input  parameters  for  a  proposed 
test  configuration.  The  control  volume  theory  will  be 
compared  with  experimental  data  in  Chap.  3  and  also  with  more 
accurate  digital  computer  model  results  in  Chap.  4. 
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CHAPTER  3 


EXPERIMENTAL  FORMULATION 

The  purpose  of  this  chapter  is  to  present  the  second  of 
three  studies,  the  experimental  formulation  of  a  cylindrical 
body  falling  through  fluid  in  a  closed-end  tube.  The  results 
of  this  experimental  effort  are  compared  to  the  results 
predicted  by  the  control  volume  formulation  and  will  be 
subsequently  compared  to  the  numerical  (differential) 
formulation. 

A  brief  example  helps  clarify  the  dynamic  scope  of  the 
experimental  problem.  The  prototype  could  consist  of  either  a 
10  foot  diameter  train  traveling  at  100  mph  (147  fps)  in  air, 
or  a  21  inch  diameter  torpedo  traveling  at  55  fps  in  water. 
In  both  cases,  Re^»Vd/v  =  9E6.  Unfortunately,  the  model  test 
could  not  achieve  this  high  Reynolds  number.  The  highest 
Reynolds  number  actually  achieved  in  testing  was  only  22700. 
As  discussed  in  the  next  two  sections,  very  few  practical 
options  are  available  to  realize  high  test  Reynolds  numbers. 


Section  3.1  Methods  of  Obtaining  High  Model  Reynolds  Numbers 
There  are  five  options  available  to  increase  Re^«=Vd/v. 
The  first  three  involve  increasing  the  terminal  velocity. 
Terminal  velocity,  V^,  from  the  theory  in  Chap.  2,  eq.  2-4, 
can  be  rewritten  as  follows; 


2A(S-I)g 


.,0.5 


a  C 


dt 


3-la 
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or  V,.  =  K  US  -  Dgl'*'*,  “ith  K  -  — - 

where  d  -  model  diameter 

2 

a  “  model  frontal  area  «  nd  /4 
A  «  model  volume 

«  equivalent  model  length  *  A/a 
§  -  specific  gravity  of  the  model 
g  «  acceleration  of  gravity 

-  drag  coefficient  at  terminal  velocity . 
Examining  eq.  3— lb  it  is  clear  that  terminal  velocity,  V^,  in 
a  drop  test  can  be  raised  by 

Option  #1:  increasing  model  specific  gravity,  S 
Option  #2;  increasing  model  equivalent  length, 

Option  #3:  decreasing  drag  coefficient, 

Option  #1  indicates  that,  all  else  being  equal,  a  model 
made  of  stainless  steel,  §=7.66,  will  have  a  terminal 
velocity  1.7  times  larger  than  a  model  made  of  aluminum, 
§=2.78.  As  discussed  later  in  more  detail,  it  was  decided  to 
use  aluminum  because  it  is  much  easier  to  machine  and  drill. 
The  choice  of  §=2.78  meant  that  option  #1  was  not  fully 
utilized  to  increase  terminal  velocity. 

Option  #2  indicates  that  models  with  a  long  effective 
length,  (i.e. , large  A  and  small  a)  will  have  the  highest 

terminal  velocity.  However,  large  values  of  model  length,  L, 
require  long  tube  lengths,  which  is  not  practical.  It  was 

also  decided  to  independently  vary  e=D/d  and  to  fix  L  such 
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that  certain  dimensionless  parameters  would  remain  constant 
as  discussed  in  Sect.  3.7,  "Dimensionless  Groups,"  later  in 
this  chapter.  Option  #2  was  therefore  not  used  to  increase 
terminal  velocity. 

Option  #3  was  implemented  by  the  use  of  a  streamlined 
model  geometry,  as  described  in  Sect.  3.6. 

Option  #4  is  to  increase  Rej>=Vd/v  through  large  model 
dimensions , i .e. , large  model  diameter,  d.  This  might  partially 
compensate  for  lower  terminal  velocities  arising  from  the 
choices  made  by  the  author  as  discussed  under  options  #1,  #2 
and  #3  above.  The  limitations  associated  with  large  model 
dimensions  were  discussed  under  option  #2. 

Option  #5  is  the  last  method  to  increase  Re^“Vd/v,  which 
is  to  decrease  the  kinematic  viscosity,  v,  of  the  test  fluid. 
This  method  would  allow  for  the  desired  relatively  small 
model  dimensions  and  low  terminal  velocities.  Given  that  the 
model  might  consist  of  a  1  in.  diameter  cylinder  traveling  at 
5  fps,  for  Re^=9E6  the  fluid  should  have  a  viscosity  on  the 
order  of  4.7E-8  ft  /s,  which  is  a  factor  of  about  225  less 
than  water.  Appendix  A  contains  several  rather  complete 
tables  of  kinematic  viscosities  for  both  gases  and  liquids 
and  discusses  methods  of  achieving  low  viscosities.  The 
conclusion  drawn  in  Appendix  A  is  that  there  is  no  simple 
practical  method  available  for  producing  such  a  low 
viscosity. 

Based  on  the  limited  options  to  significantly  raise 
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Re^,  given  the  choice  of  small  streamlined  aluminum  models 
traveling  at  low  terminal  velocities,  it  was  decided  to  use  a 
simple  drop  test  consisting  of  a  wire-guided  aluminum  model 

with  water  as  the  working  fluid. 

The  only  other  readily  available  experimental  option  open 

to  the  author  was  to  consider  wind  tunnel  testing.  For 
reasons  discussed  in  the  next  section  wind  tunnel  testing  was 
not  used  during  this  study. 


section  3.2  Feasibility  of  Testing  in  a  Wind  TunneJ. 

A  brief  calculation  shows  that  in  a  wind  tunnel (with  V«100 
ft/s,  d-3.0"  and  v=1.65E-4  ft^/s)  Re^.Vd/v=l . 5E5, which  is 

about  6.7  times  larger  than  the  maximum  Re ^=22700  achieved  in 
the  drop  test.  While  this  larger  Reynolds  number  is  obviously 
desirable,  there  are  two  main  reasons  why  wind  tunnel  testing 
was  not  done. 

Reason  #1.  Violation  of  Boundary  Conditions 
For  kinematic  similitude  to  exist, one  must  arrange  for  a 
moving  wall  in  axisymmetric  flow,  which  is  essentially 
impossible.  Therefore,  wind  tunnel  data  must  be  taken  with 
stationary  tube  walls.  These  data  will  be  reasonably  good 
when  e-D/d  is  large,  in  which  case  the  streamline  pattern 
around  the  model  is  not  adversely  affected  by  the  boundary 
layer  growth  on  the  tube  walls  due  to  the  bulk  flow  in  the 
tunnel.  In  the  prototype  there  is  no  bulk  flow  through  the 
tube.  For  the  small  values  of  3=1.14,  1.33  and  1.6  of 
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interest  to  this  study,  streamline  displacement  is  sufficient 
to  alter  drag. 

Reason  #2.  Distorted  Approach  Velocity  Profile 

In  the  prototype,  the  boundary  layer  on  the  tube  wall  due 
to  the  moving  body  is  confined  essentially  to  the  annular  gap 
between  the  body  and  the  wall.  The  fluid  is  at  rest 
everywhere  in  the  field  except  in  the  near  field  surrounding 
the  body.  However,  in  a  wind  tunnel,  the  bulk  flow  through  the 
tunnel  causes  a  boundary  layer  to  build  along  the  entire 
length  of  the  tube  wall.  At  the  point  where  the  bulk  flow  in 
the  wind  tunnel  accelerates  into  the  annular  gap,  the 
presence  of  the  tube  boundary  layer  distorts  the  velocity 
profile  of  the  flow  past  the  model.  This  distortion  will  be 
less  objectionable  for  large  P's  but  for  the  small  clearances 
of  interest  to  this  study  the  profile  distortion  will  produce 
unreliable  drag  results. 

The  rest  of  this  section  discusses  the  tube  boundary  layer 
in  a  wind  or  water  tunnel.  Estimates  of  the  streamline 
displacement  and  tube  boundary  layer  height  are  made  and 
compared  to  the  gap  height,  h=(D-d)/2.  A  table  is  presented 

c: 

summarizing  sample  calculations. 

Fluid  is  displaced  away  from  the  tube  wall  by  a  distance, 
6^  ,  called  the  displacement  thickness.  Estimate  on  the 

tube^  wall  using  flat  plate  turbulent  flow  relations  from 
White  [39],  figure  5  ; 


1.  subscript:  T  =  Tube 
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3-2a 


« 


z 


.  T 

where  - 

z 


8 


/ 


z 


0.16 


3-2b 


for  V  -  free  stream  velocity  entering  the  tube 

6^  -  boundary  layer  thickness  on  the  tube  wall,  defined 
to  be  the  radial  distance  from  the  wall  where  the 


fluid  velocity  is  equal  to  99%  of  V 
5^*  .  displacement  thickness  on  the  tube  wall  -  S,j,/8 
z  -  axial  distance  along  the  tube  measured  from  the 
entrance  of  the  tube 

V  »  kinematic  viscosity  of  the  fluid  medium. 
Multiply  eq.  3-2b  by  z/D; 


★ 

=  0.02 

'z' 

6/7 

'VD" 

D 

.D. 

.V  . 

-1/7 


3-3a 


VD  ved 

where  Re_  *  —  “  “  P  ’ 

V  V 

Multiply  eq.  3-2a  by  z/D: 


D  D 

Examine  the  gap  clearance,  h,  where  h=(D-d)/2: 

h  e-1 


D 


Then , 


h 


2P 

h/D 


eq.  3-3a 


eq.  3-4 


3-3b 


3-3c 


3-4 


3-5 
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3-6 


Multiply  eq.  3-3c  by  D/h  to  get: 


h  h 


Representative  calculations  from  eq.  3-3  through  eq.  3-6  for 
a  wind  tunnel  and  water  tunnel  are  summarized  in  Table  1  . 


Table  1.  Displacement  Thickness 
Wind  Tunnel  and  Water  Tunnel 


Variable 

Equation 

Ai  r 

Water 

d 

3.0  in. 

3.0  in 

e 

1.6 

1.6 

D 

4.8  in. 

4.8  in 

Z/D 

5. 

5. 

V 

1.65E-4  ft^/s 

1.07E-5  ft  Vs 

V 

100.  fps 

20.  fps 

ReT 

3-3b 

2.4E5 

7.5E5 

3-3a 

1.35% 

1.15% 

3-3c 

10.8% 

9.2% 

h/D 

3-4 

18.8% 

18.8% 

3-5 

7.2% 

6.1% 

3-6 

57.6% 

49.1% 

Thus, 

for  the  values 

of  d,e,D,z/D  and 

V  of  interest,  the 

water  tunnel  tube  boundary  layer  will 

occupy  49%  of  the 

annular 

gap  and  the  wind  tunnel  57%. 

The  inlet 

profile 

distortion  is 

apparently  quite 

extensive.  As  8  is 

further 

decreased  from 

1.6  down  to  1.14, 

the  inlet  velocity 

profile  engulfs  the  gap  and  there  is  total  distortion. 
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The  preceding  considerations  highlight  the  extremely 
difficult  task  of  obtaining  model  data  at  prototype  Reynolds 
numbers.  The  next  section  covers  the  actual  model  test 
consisting  of  a  drop  test,  using  water,  with  a  short  tube  and 
three  small  wire-guided  models. 

Section  3.3  Drop  Test  Apparatus 

The  model  started  from  rest  and  did  a  free-fall  inside  a 
water  filled  glass  tube,  figure  7  .  Concentricity  was 
ensured  by  the  use  of  a  guidewire  passing  down  the  centerline 
of  the  tube  and  threaded  along  the  model  centerline.  The  tube 
was  made  by  Corning  Glass  Works  and  had  a  1"  inside  diameter 
and  was  6  ft.  long  (D-1.0",  . 0  ft.). 

The  frame  supporting  the  tube  was  built  so  that  the  tube 
could  be  made  absolutely  vertical  to  ensure  that  the  body 
fell  straight  down  to  minimize  any  normal  force  (hence 
sliding  friction)  between  the  guidewire  and  the  model.  The 
frame  was  set  on  a  pin  at  its  base  and  then  adjusted 
f ront-to-back  by  a  set  of  turnbuckles.  It  was  also  adjusted 
right-to-left  by  a  second  set  of  turnbuckles.  It  was 
determined  to  be  vertical  by  means  of  a  carpenter's  level 
held  to  the  side  of  the  tube. 

At  the  bottom  of  the  tube  a  hard  rubber  stopper  with  a 
small  central  hole  for  the  guidewire  was  sealed  into  place  by 
means  of  room  temperature  vulcanizing  sealant  (RTV).  The  top 
of  the  tube  also  had  a  similar  hard  rubber  stopper  with  a 
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small  central  hole,  which  was  firmly  seated  into  the  tube  but 
was  not  sealed  with  RTV  because  it  needed  to  be  removed 
periodically  to  change  models.  The  guidewire  was  20  lb 
monofilament  line  that  was  stretched  taut  down  the 
centerline  of  the  tube  by  passing  through  the  central  holes 
of  the  hard  rubber  stoppers  on  either  end  of  the  tube. 
Another  very  small  hole  was  made  off  the  centerline  in  the 
top  rubber  stopper  to  allow  a  very  small  diameter  (eO.lmm) 
piece  of  monofilament  line  to  pass  through  so  that  the  model 
could  be  pulled  back  up  the  tube  for  the  next  run.  This 
retrieval  line  was  held  by  a  solenoid  operated  plunger  that 
was  actuated  at  the  beginning  of  each  run  to  release  the 
model.  The  solenoid  mechanism  and  its  power  supply  were  built 
by  the  author  out  of  standard  parts  available  at  retail 
electronics  stores. 

Section  3.4  Determination  of  Terminal  Velocity 

Various  techniques  were  evaluated  for  determining  terminal 
velocity,  since  it  is  such  a  critical  parameter.  In  all  cases 
the  position  of  the  model  as  a  function  of  time  was  recorded. 
Velocity  was  inferred  from  accurate  knowledge  of  distance  and 
time.  Using  the  constant-C^  model  from  Chap.  2,  it  was 
determined  that  the  upper  limit  of  terminal  velocities  was 
going  to  be  about  5  fps.  All  position  vs.  time  data  presented 
in  this  report  were  taken  using  a  high  speed  movie  camera  set 
at  100  frames/sec.  The  test  set-up  is  schematically  shown  in 
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figure  7  . 

The  high  speed  movie  camera  had  a  built-in  timing  light 
and  gave  very  repeatable  position  vs.  time  results,  as 
presented  in  Sect.  3.10.  It  was  attempted  to  get  model 
location  using  an  aluminum  yardstick  mounted  to  the  wooden 
frame  behind  the  glass  tube  containing  the  model.  This  method 
ended  up  providing  only  a  rough  check  since  the  fine  lines  on 
the  yardstick  turned  out  to  be  unreadable  on  several  of  the 
films.  Precise  location  of  the  model  was  determined  by  the 
technique  described  in  Sect.  3.9.  Elapsed  time  was  precisely 
measured  using  a  strobe  light  as  a  backup  to  the  internal 
camera  timing  light. 

Section  3.5  Determination  of  High  Speed  Movie  Camera  Settings 

Reference  is  made  to  the  case  discussed  in  Sect.  2.11 
"Application  to  Design  of  the  Experiment."  Equation  2-4  is 
used  to  calculate  time,  T«0.641  sec,  required  to  reach  98%  of 
terminal  velocity.  This  is  the  time  required  to  fall  from 
rest  to  23"  of  displacement.  Also,  solve  =  53  in./s  from  eq. 
2-4,  for  use  in  determining  camera  speed.  Through  trial  and 
error  it  was  determined  the  model  should  move  about  0.5 
inches  from  one  frame  to  the  next  frame  to  keep  blurring  to 
an  acceptable  minimum.  Thus  the  camera  was  set  at  100  frames 
per  second.  This  setting  happens  to  produce  a  convenient  0.01 
sec/frame.  The  camera  required  1  second  to  start  up  and 
reach  a  steady  state  film  speed.  Since  T=0.641  sec,  another 
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second  was  allowed  to  run  the  test,  and  one  final  second 
provided  for  the  camera  to  coast  to  a  stop.  Thus,  total  camera 
time  was  3  seconds  or  about  300  frames  per  run.  At  about  12 
frames  per  foot  of  film  length,  this  requirement  translated 
into  25  feet  of  film  per  run.  It  also  hints  at  the  physical 
difficulty  experienced  manipulating  large  quantities  of  film 
to  obtain  position  and  time  data. 

Section  3.6  Model  Construction  and  Dimensions 

When  building  the  models  there  was  a  problem  drilling  a 
very  small  diameter  (0.02")  hole  down  the  length  of  the  model 
to  allow  the  guidewire  to  pass  through,  because  small 
diameter  drills  have  a  very  strong  tendency  to  "wander"  off 
the  centerline  and  also  to  break.  This  problem  was  resolved 
by  drilling  a  larger  hole,  dj^,  down  the  centerline  from  each 
end,  meeting  in  the  middle.  The  larger  drill  "wandered"  far 
less  and  didn't  break  like  the  small  drills.  Then  at  each  end 
of  the  model  a  small  bead  was  welded  onto  the  model  covering 
the  1/4"  hole.  The  weld  bead  was  then  drilled  to  0.02"  to 
allow  for  the  guidewire  to  pass  through  the  centerline  of  the 
model.  It  was  felt  that  this  arrangement  1)  created  the 
least  friction  between  the  guidewire  and  the  model  since 
there  was  essentially  only  point  contact  at  the  ends  vs.  line 
contact  through  the  entire  model,  2)  insured  that  the  model 
was  guided  on  its  centerline,  3)  was  practical  to  implement. 

An  additional  0.017"  diameter  hole  was  drilled  in  each 
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model  to  accommodate  the  retrieval  line.  To  minimize  any 
adverse  effects  on  the  flow  field,  this  hole  was  drilled 
perpendicular  to  the  centerline  about  3/8"  from  the  apex  of 
the  cone,  on  the  tail  of  the  model. 

Summarizing,  figure  8  ,  three  bodies  0=1.6,  3=1.33  and 
3-1.14)  were  tested  to  get  distance  vs.  time  data  while 
falling  inside  a  1.000"  tube,  figure  7  .  They  were  called 
625,  750  and  875,  which  corresponds  to  outside  diameters  of 
5/8,  3/4  and  7/8  inches,  respectively.  All  the  bodies  had  a 
hemispherical  nose,  cylindrical  mid— body,  and  a  conical  tail 
(also  referred  to  as  the  nose,  straight  and  tail  sections 
respectively).  Table  2  lists  relevant  model  geometries  and 
properties,  which  are  sketched  in  figure  8  .  Note  that  A  is 
the  volume  of  water  the  model  displaced  in  a  ready-to-run 
condition  and  that  terminal  velocities  are  the  experimentally 
obtained  values. 


50 


Table  2. 

Model 

2 

Specifications 

Body  625 

Body  750 

Body  875 

Model  Parameters 

Diameter ,  d  (  in. ) 

0.625 

0.750 

0.875 

Center  hole,  dj^  (in.) 

0.125 

0.125 

0.125 

(in.) 

L_  (in.) 

0.313 

0.375 

0.438 

5.433 

3.453 

1.588 

s 

L.  (in.) 

1.38 

1.375 

1.485 

L»L  +L  +L.  (in.) 

7.125 

5.203 

3.510 

n  s  t  2 

Frontal  area,  a  (in.  ) 

0.307 

0.442 

0.601 

Displaced  Vol .  ,  A  (in.^) 

1.64 

1.79 

1.37 

Lg  -  A/a  (in.) 

5.33 

4.05 

2.27 

Dry  weight,  W  (lb) 

0.1640 

0.1796 

0.1377 

Buoyant  force,  F^^  (lb) 

0.0591 

0.0646 

0.0494 

Net  wt,  -  W-Fj^  (lb) 

0.1050 

0.1150 

0.0884 

Term.  Vel.,  V^,  (in/s) 

55.1 

30.5 

9.0 

Primary  Dimensionless  Groups  (evaluated  at  V^.) 

»  A/D^ 

1.64 

1.79 

1.37 

II3  -  9D/V^2 

0.123 

0.413 

4.88 

=  Re^  -  V^.D/v 

3.63E4 

2.E4 

5.77E3 

n  »  e  »  D/d 

1.600 

1.333 

1.143 

Hg  -  5 

2.78 

2.78 

2.78 

Secondary  Dimensionless 

Groups 

(evaluated  at  ) 

”l  ■ 

0.358 

1.32 

11.9 

"7  “  Rejj-v^d/v-n^/ng 

2.27E4 

1.5E3 

5.E3 

Hq  -  Lg/d  -  (n5^n2)4/n 

8.53 

5.40 

2.60 

«9  - 

2.34 

5.97 

39.6 

A  discussion  of  the 

dimensionless  groups  0^ 

through  Ilg 

presented  in  the  next  section. 


2.  Fluid  is  water,  p  «  1.938  slugs/ft^,  v=1.07E-5  ft^/sec 
fi  “  2.074E-5  slug/ft-sec 


is 
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Section  3»7  Dimensionless  Groups 

The  problem  of  a  model  falling  in  a  tube  can  be 

functionally  described  by  f ( F^, S , A,g, p,^,D,d,V)-0,  where 

2 

-  drag  force  on  the  model  (ML/t  ) 

S  -  specific  gravity  of  the  model  (1) 

A  “  displaced  volume  of  the  model  (L^) 

2 

g  “  acceleration  of  gravity  (L/t  ) 
p  -  density  of  fluid  (M/L^) 

-  dynamic  viscosity  of  fluid  (M/Lt) 

D  »  tube  diameter  (L) 

d  -  model  diameter  (L) 

V  -  model  velocity  (L/t). 

Using  p,V,D  the  following  dimensionless  groups  appear; 

"l  “  F^/(pV^D^) 

112  *  A/D^ 

"3  “  gD/v^ 

-  Re^  -  pVD/// 

Hg  =  e  -  D/d 

h  -  s- 

Functionally,  the  problem  may  be  expressed  as  follows: 

It  is  clear  that  112,  ^6  fixed  through  the  choice  of 

a  model  and  a  tube.  For  this  study, 

112  "  1-65  for  all  cases.  This  value  is  of  the  greatest 

interest  to  the  Navy,  corresponding  to  a  21"  diameter 
torpedo  in  a  27"  diameter  tube.  The  actual 
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experimental  value  varied  +15%  due  to  the  problem  of 
drilling  the  central  guidewire  hole  per  Sect.  3.6. 

D  “  1.00"  fixed  tube  diameter. 

S  -  2.78,  aluminum  models. 

d  -  0.625",  0.750",  0.875"  fixed  model  sizes, 

corresponding  to  8*1 . 6 , 1 . 33 , 1 . 14  defined  by  specific 
torpedo-tube  combinations  of  interest  to  the  Navy. 
Examining  it  appears  that  the  dimensional 

problem  is  to  determine  drag  force,  F^,  as  a  function  of 
velocity,  V,  for  a  particular  model  (S,d,A)  in  a  fixed  tube 
(D)  and  predetermined  fluid 

It  is  possible  to  recast  the  dimensionless  problem  into 
more  familiar  terms. 

Let  -  Re^  «  Vd/v  -  , 

and  Dg  =  L^/d  «=  A/ad  =  {ng^n2)4/ii. 

At  terminal  velocity,  V^,  it  is  known  that  i.e., 

drag  force  equals  the  net  weight  of  the  model,  where 

Wjj  “  gp(S-l)A  -  n2n3(ng-l)pV^D^, 

or  -  W^/(pV^^D^)  -  n2n3(ng-l)  ,  at  only. 

It  is  customary  to  define  a  drag  coefficient,  in  terms  of 

2 

model  velocity,  V,  and  model  frontal  area,  a*d  r/4: 

flg  -  Cd  -  2F^/(pV^a)  ,  at  all  times, 
and  Hg  -  (nj^n5^)8/n  ,  at  terminal  velocity  only. 

The  dimensionless  problem  then  becomes  finding  the  function, 
h,  where 

«  h(Re^,8,Lg/d,fluid), 
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or 


ng“h(n2»n2/n^jng,ng). 

Experimentally,  only  one  point  of  this  function  was 
obtained  for  each  model.  This  point  corresponded  to  the 
terminal  velocity  where  the  exact  value  of  the  drag  force, 
F^,  is  known.  As  discussed  further  in  Chap.  4,  the  function 
h  was  predicted  using  the  constant— model  and  compared  to 
predictions  made  by  the  more  accurate  numerical  model. 

Section  3.8  Test  Procedure 

A  typical  run  proceeded  as  follows: 

1.  Initialize  the  model.  Using  the  retrieval  line  the 
model  would  be  gently  pulled  back  up  the  tube  to  a  preset 
position,  at  the  top  of  the  camera  window.  At  this  point  the 
retrieval  line  would  be  placed  under  the  solenoid  operated 
release  mechanism. 

2.  Lights.  Proper  light  levels  were  critical  to  prevent 
under  or  overexposure  of  the  film.  Due  to  the  short  time  each 
frame  was  exposed,  high  powered  lights  had  to  be  used. 

3.  Camera.  The  high  speed  movie  camera  required  about  1 
second  to  come  up  to  speed,  so  at  this  point  the  camera 
switch  was  manually  activated. 

4.  Action.  The  model  was  released  by  activating  the 
solenoid  that  held  the  retrieval  line.  At  the  same  time  a 
strobe  light  started  up  to  mark  the  beginning  of  the  run  on 
the  film. 

5.  End  of  run.  Camera  and  lights  were  shut  off. 
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6.  Record  run  number  and  make  notes  of  camera  speed,  body 
number,  etc. 

The  camera  had  to  be  focused  very  precisely  or  else  the 
nose  of  the  model  was  in  focus  while  the  side  of  the  model 
closest  to  the  camera  was  blurry.  Focus  was  an  important 
consideration  when  reducing  the  data  as  described  below. 

Section  3.9  Reducing  the  Data 

Data  were  obtained  by  examining  the  film  using  a  microfiche 
machine.  This  allowed  the  developed  film  to  be  viewed  one 
frame  at  a  time  on  an  enlarged  screen.  The  nose  of  the  model 
had  to  be  in  sharp  focus  and  the  level  of  film  exposure  was 
critical.  A  measurement  was  taken  of  the  body  length  on  the 
screen,  which  was  compared  to  the  precisely  known  actual 
model  length,  so  that  a  scale  could  be  made.  Then  for  each 
frame  the  distance  the  nose  moved  from  the  previous  frame, 
called  SS,  was  measured  to  within  0.01  inch  and  was  recorded. 
Each  frame  also  had  a  timing  mark  on  the  border  recorded  by 
the  internal  strobe  light  in  the  camera.  Elapsed  time  from 
the  previous  frame,  called  8^,  was  interpolated  and  recorded. 
This  process  continued  frame-by-frame  until  the  body 
disappeared  off  the  screen.  To  ensure  accuracy,  the  film  was 
rewound  to  the  beginning  of  the  run  and  the  whole  process  was 
repeated  frame-by-frame.  At  the  conclusion  of  the  second 
reading  the  results  were  compared  and  differences  were 
resolved  by  checking  any  frames  in  question.  It  was  felt  that 
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this  procedure  minimized  self-error. 


Section  3.10  Presentation  of  Test  Results 

After  the  data  for  each  run  had  been  tabulated  into 
columns  of  8S  and  8T,  as  described  above, they  were  entered  in 
data  files  on  the  VAX  8600  super-minicomputer  to  facilitate 
graphing  and  presentation.  It  was  found  that  the  dataware 
vgcy  repeatable  for  all  three  bodies.  Tabulated  displacement 
vs.  elapsed  time  results  for  all  three  models  are  contained 
in  Appendix  A.  The  625  model  results  are  contained  in  Table 
A-4  .  Table  A-5  is  a  tabulation  of  results  for  the  750 
model  with  Table  A-6  tabulates  results  for  the  875  model. 

Table  3  presents  a  summary  of  time  vs.  average 
displacement  for  each  of  the  three  bodies,  625,  750  and  875. 
Asterisks  are  given  for  displacement  when  the  model  is  at 
terminal  velocity.  For  instance,  itwas  estimated  that  the  875 
model  was  at  terminal  velocity  after  it  had  traveled  2.05". 
However,  the  750  model  required  10.67  inches  and  the  625  a 
full  19.63  inches  to  reach  terminal  velocity.  These  results 
are  compared  to  the  constant— model  composed  of  eq.  2—5  and 
eq.  2-20,  in  Sect.  3.11. 

There  had  been  some  concern  that  the  model  would  fall  at 
slightly  different  rates  due  to  very  small  changes  in  angle 
of  attack,  0,  as  it  went  down  the  guidewire,  but  there  was  no 
discernible  hint  of  this  in  the  data.  It  is  probably  a 
secondary  effect. 
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Table  3.  Averaged  Displacement  of  Models  vs- 


Time 

sec . 

Model 

625 

Z 

avg 
in . 

Model 

750 

^avg 
in . 

Model 

875 

Z 

avg 
in . 

0.00 

0.00 

0.00 

0.00 

0.01 

0.02 

0.01 

0.01 

0.02 

0.06 

0.03 

0.03 

0.03 

0.12 

0.06 

0.07 

0.04 

0.19 

0.11 

0.13 

0.05 

0.29 

0.19 

0.19 

0.06 

0.41 

0.29 

0.25 

0.07 

0.56 

0.39 

0.31 

0.08 

0.74 

0.51 

0.38 

0.09 

0.91 

0.62 

0.46 

0.10 

1.09 

0.74 

0.56 

0.11 

1.30 

0.89 

0.64 

0.12 

1.55 

1.06 

0.72 

0.13 

1.77 

1.25 

0.80 

0.14 

2.03 

1.45 

0.89 

0.15 

2.29 

1.63 

0.97 

0.16 

2.59 

1.83 

1.06 

0.17 

2.90 

2.08 

1.14 

0.18 

3.19 

2.31 

1.23 

0.19 

3.52 

2.53 

1.31 

0.20 

3.90 

2.78 

1.40 

0.21 

4.27 

3.05 

1.50 

0.22 

4.67 

3.28 

1.59 

0.23 

5.07 

3.55 

1.68 

0.24 

5.47 

3.78 

1.77 

0.25 

5.87 

4.07 

1.87 

0.26 

6.29 

4.33 

1.96 

0.27 

6.68 

4.59 

2.05 

0.28 

7.14 

4.86 

*** 

Time 
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Table  3.  continued 


Model 

Model 

Model 

625 

750 

875 

Time 

sec . 

z 

avg 
in . 

^avg 
in . 

^avg 
in . 

0.29 

7.58 

5.15 

★  ★  ★ 

0.30 

8.03 

5.44 

*  ★  * 

0.31 

8.52 

5.73 

*  ★  ★ 

0.32 

9.05 

6.00 

★  ★  ★ 

0.33 

9.51 

6.27 

★  ★  ★ 

0.34 

9.97 

6.56 

★  ★  ★ 

0.35 

10.48 

6.86 

*  *  ★ 

0.36 

10.98 

7.15 

*  *  * 

0.37 

11.47 

7.44 

★  *  * 

0.38 

11.99 

7.73 

*  ★  * 

0.39 

12.49 

8.02 

0.40 

13.01 

8.31 

★  *  ★ 

0.41 

13.53 

8.60 

ic  if  it 

0.42 

14.06 

8.89 

it  it  it 

0.43 

14.59 

9.19 

it  it  it 

0.44 

15.13 

9.50 

it  it  it 

0.45 

15.69 

9.79 

*  *  * 

0.46 

16.25 

10.08 

*** 

0.47 

16.82 

10.37 

*  ★  * 

0.48 

17.40 

10.67 

★  ★  ★ 

0.49 

17.97 

*** 

it  it  ^ 

0.50 

18.54 

*** 

it  it  it 

0.51 

19.11 

*** 

it  it  it 

0.52 

19.63 

*** 

it  it  it 

0.53 

**  * 

it  ir  it 

it  it  it 

The  displacement  vs.  time  data  for  the  625  model  have  been 
plotted  in  figure  9  .  It  is  seen  from  figure  9  that  the 
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data  correlate  very  closely  to  the  averaged  displacement, 

Z  The  error  bands  on  figure  9  embrace  the  maximum  and 

avg’ 

minimum  displacement  values  for  any  individual  run.  The  data 
are  totally  contained  within  a  ±4%  band  around  with  a 

standard  deviation  of  ±3%.  Figure  9  also  compares  the 
constant-Cjj  prediction  made  from  eq.  2-5  with  Z^^g.  Eq.  2-5 
is  consistently  low,  but  at  the  worst  is  only  off  by  6%.  It 
is  thought  that  the  constant-C^  model  under-predicts 
displacement  because  the  use  of  over-predicts  drag  in 

the  early  stages  of  the  run,  when  the  velocity  of  the  model 
is  very  small. 

A  variable-C^  model  was  discussed  in  Sect.  2.12  where 
was  allowed  to  vary  with  velocity,  figure  18  ,  as  opposed  to 
being  fixed  at  the  terminal  velocity  value,  results 

of  the  variable— model  as  computed  by  the  Runge— Kutta  based 
computer  program  in  Appendix  E.l  are  also  shown  in  figure 
9  for  comparison  to  data  and  the  constant— model  from  eq. 
2-5.  It  is  seen  from  figure  9  that  the  variable-C^  model 
consistently  under-predicts  the  data  and  eq.  2—5,  being  as 
much  as  8%  lower  than  at  terminal  velocity.  A  possible 

civ  g 

explanation  of  this  result  is  that  the  variable-C^  model 
assumes  turbulent  flow  for  the  entire  run,  i.e.,  from  release 
through  terminal  velocity.  The  consequence  of  this  assumption 
is  that  drag  is  over-predicted  during  the  initial 
displacement,  causing  the  variable-C^j  model  to  under-predict 
displacement.  This  under-prediction  is  slightly  worse  than 
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the  constant-C^  model  because  the  variable-C^  model  predicts 

Cj>Cj  up  to  the  terminal  velocity.  It  is  concluded  that  the 
d  dt 

constant-C^  model  is  acceptable  and  preferred  over  the 

variable-Cj  for  routine  engineering  calculations  in  the 
d 

startup  region. 

Figure  10  graphically  displays  test  results  for  the  750 
model.  The  maximum  and  minimum  data  are  enti rely  contained  in 
an  envelope  +4%  around  with  a  standard  deviation  of  +3%. 

Figure  11  presents  the  test  results  for  the  875  model. 
The  maximum  and  minimum  displacement  for  any  particular  run 
are  found  to  be  contained  in  a  +5%  band  with  a  standard 
deviation  of  +4%  with  reference  to 

In  Chap.  4  the  constant-C^  model  will  be  compared  to  the 
numerical  model  over  a  wide  range  of  Reynolds  numbers  for 
P-1.6,  1.33  and  1.14. 

Section  3.11  Curve-Fit  Model 

Although  the  constant-C^  model  adequately  predicted 
displacement  of  the  model  as  discussed  in  Sect.  3.10,  it  is 
nice  to  have  a  little  more  accuracy  during  the  acceleration 
phase  of  travel  in  order  to  more  accurately  predict  velocity 
and  displacement.  It  is  proposed  to  use  a  least— squares— curve 
fitting  process  to  determine  the  constants  a^  and  b^  in  the 
constant— C^  model,  eq.  2—4  and  eq.  2—5.  The  coefficients  a^ 
and  b  will  be  called  a„^  and  b„^  respectively. 

Table  4  contains  values  of  a^,  b^,  and  C^^  from  the 
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constant— model  of  eq.  2-4,  eq.  2-5  and  eq.  2—20.  Table  4 
also  lists  the  least-squares-curve-f i t  a^^  and  b^^.  It  is 
noted  that  the  curve-fit  for  a^^  and  b^^  is  very  stiff,i.e.,  a 
fairly  wide  range  of  values  for  a^^  and  b^^  will  fit  the 
data.  Values  for  Re,  and  Re„  come  from  eq.  2-14  at  steady 
state  velocity.  Values  for  given  in  eq.  3-9  later  in 

this  section,  are  given  for  convenience  to  provide  a 
comparison  to 


Table  4.  Comparison  of  Measured  and  Predicted  Values 


Model 

^oe 

^o 

inches 

^e 

sec”^ 

^e 

'^t 

in/s 

^dte 

^dt 

Rea 

625  Model 

measured 

13.67 

4.03 

55.1 

2.34 

22700. 

22400. 

predicted 

15.61 

3.59 

56.0 

2.26 

750  Model 

measured 

6.09 

5.01 

30.5 

5.94 

14800. 

11300. 

predicted 

5.55 

5.52 

30.6 

5.90 

875  Model 

measured 

0.715 

12.6 

9.0 

40.1 

5100. 

3100. 

predicted 

0.686 

12.9 

8.9 

41.5 

Figure 

9  compares  to 

the  curve-fit  model  for  the 

625  model 

.  The 

curve-fit 

model 

varies 

only  +1% 

^avg* 

Figure  10  presents  and  the  curve-fit  for  the  750 


61 


model,  with  a  +1%  deviation.  Figure  11  shows  the  curve  fit 
model  varying  +2%  from  for  the  875  model.  Once  the 

curve— fit  model  has  been  verified, it  is  possible  to  make  very 
accurate  predictions  of  displacement,  velocity  and 
acceleration  during  the  acceleration  phase,  using  a^^  and  b^^ 
in  eq.  2-4  and  2-5.  This  has  been  done  for  the  625  model  and 
the  results  are  given  in  figure  15  .  For  example,  referring 
to  figure  15  ,  when  t=0.3  s  the  acceleration  of  the  625 

model  is  about  39in./s^,  velocity  is  roughly  46  in./s  and  the 
model  has  moved  a  total  of  approximately  8  inches  from  the 
time  it  was  released.  At  t=0.6  s,  the  625  model  is  at  54in./s 
and  has  traveled  23  in.  This  result  was  predicted  in  Sect. 
2.11  using  the  constant-C^  model. 

Figure  15  is  also  useful  for  determining  when  the  625 

model  is  at  terminal  velocity.  It  is  noted,  from  eq.  2-3, 

that  at  the  moment  the  model  is  released,  acceleration  is 

always  equal  to  Wj^/(m+mj^).  In  the  next  instant  after  release, 

when  the  model  has  a  velocity  greater  than  zero,  the 

acceleration  begins  to  fall  due  to  drag.  Thus,  it  is  possible 

to  establish  a  criterion  for  terminal  velocity  defined  as  the 

model  velocity  when  acceleration  goes  below  O.Olg  or  about  4 

in./s^ .  Looking  at  figure  15  it  is  seen  that  the 

.  2 

acceleration  of  the  625  model  has  fallen  under  4  in./s  at 
t-0.6  sec.;  similarly  for  the  750  model,  figure  16  predicts 
terminal  velocity  of  about  30.5  in./s  occurring  at  t=0.46  sec 
when  displacement  is  about  11  inches.  The  875  model  results 
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are  given  in  figure  17  where  the  terminal  velocity  of  9 
in./s  occurs  when  t=0.3  sec.  and  displacement  is  a  little  over 
2  inches. 

Once  data  are  available  it  is  possible  to  get  an 
experimentally  determined  hydrodynamic  mass  and  drag 
coefficient  from  the  curve-fit  constants  a^^  and  b^^ .  Solve 
eq.  2-4a  for  hydrodynamic  mass  and  eq.  2-4b  for  terminal  drag 
coefficient,  called  mj^^  and  in  terms  of  a^^  and  b^^ . 


Using  eq.  2-4; 


te 


a  -b 
oe  oe  ' 


3-7 


w. 


m 


n 


he 


^t^oe 


w 


3-8 


2W 


n 


'dte  2 

paVt^ 


3-9 


2  2 

From  Appendix  D,  m^  =  m^  ( 3  -  1)  ,  3-10 

where  m  is  the  mass  of  the  fluid  in  the  annular  gap  created 

g 

by  the  cylindrical  section  between  the  model  and  the  tube 
wall.  Comparing  eq.  3-9  with  eq.  2-3  indicates  that  the  only 
difference  between  and  ^dt  the  choice  of  terminal 

velocity  to  be  used,  i.e.,  or  respectively.  Table  5 
compares  values  for  mj^  found  using  simple  inviscid  theory, 
from  eq.  2-4c,  with  values  of  mj^^  from  eq.  3-8.  Because  of 
the  stiffness  of  the  equations  resulting  in  a  wide  range  of 
acceptable  values  of  a^^  and  b^^  causing  a  wide  range  in  m^^^. 


oe 


it  is  only  possible  to  conclude  qualitatively  that  as  3 
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decreases  approaches  Values  for  nig,  the  mass  of  fluid 
in  the  gap,  found  using  eq.  3-10,  are  given  for  comparison 
purposes . 


Table  5-  Comparison  of  Inviscid  Hydrodynamic  Mass  with 
Experimentally  Determined  Values 


Equation 

2-4c 

00 

1 

rn 

3-10 

3-9 

m 

®h 

"^he 

"‘g 

“he/” 

Model 

slugs 

slugs 

slugs 

slugs 

% 

% 

^dte 

625 

5.09E-3 

1.18E-3 

4.93E-4 

2.87E-3 

23.2 

9.7 

2.34 

750 

5.58E-3 

2.58E-3 

3.45E-3 

1.56E-3 

46.2 

61.9 

5.94 

875 

4 j28E-3 

5.00E-3 

5.08E-3 

4.70E-4 

117. 

119. 

40.1 

Since  p  is 

minimum 

for  the 

875  model 

and  maximum 

for  the 

625 

model.  Table  5 

clearly  shows  that 

as  the 

gap  decreases 

in  annular  height  ( 3»*1 )  the  hydrodynamic  mass  of  the  model 
increases.  This  is  due  to  the  requirement  that  the  fluid 
be  accelerated  to  a  higher  gap  velocity,  Vg,  where  from 
continuity  (A-a)Vg  »  aV,  which  gives 


V 


which  in  turn  requires  a  greater  pressure  drop. 
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Section  3.12  Experimental  Accuracy 


Details  of  the  calculations  used  to  determine  the  level  of 
experimental  accuracy  are  provided  in  Appendix  C.  The  results 
are  summarized  in  Table  6  for  convenience: 

Table  6.  RMS  Experimental  Errors 


Variable 

RMS  error 

s,  z 

1.4% 

''te 

3.4% 

^dte 

7.2% 

Red 

3.4% 

e 

0.5% 

d 

0.2% 

a 

0.4% 

D 

0.5% 

A 

1.0% 

w 

0.8% 

Wn 

2.4% 

P 

0.5% 

V 

0.5% 

A 

2.2% 

m 

0.8% 

% 

2.5% 

®o 

8.6% 

cr 

o 

11.2% 

®oe '  ^oe 

25.0% 

%e 

25.4% 
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It  is  recommended  that  future  investigators  attempt  to 

determine  terminal  velocity,  as  accurately  as  possible, 

since  errors  in  V.  cause  considerable  problems  estimating 

te 

^dte* 

Because  the  model  tends  to  approach  asymptotically, 

there  is  always  a  very  slight  acceleration  at  the  end  of  any 
run.  The  present  study  was  carried  out  using  extreme  care  to 
minimize  the  error  in  This  was  done  by  using  great  care 

measuring  AS  over  the  greatest  possible  time,  At,  using  data 
from  the  last  several  frames  of  a  run. 

Summary 

This  chapter  has  presented  the  experimental  formulation  of 
a  cylindrical  model  falling  concentrically  in  a  vertical 
tube.  It  is  shown  that  model  Reynolds  numbers  (=2E4)  were 
more  than  two  orders  of  magnitude  less  than  the  prototype 
(=9E6).  Dimensional  analysis  showed  that  the  fundamental 
problem  was  to  determine  the  variation  of  with  model 

Reynolds  number.  Test  results  from  the  drop  test  of  three 
models,  the  625,  750  and  875,  have  been  presented.  The  test 
results  show  very  repeatable  displacement  data  vs.  time  for 
each  model  with  a  maximum  +5%  spread  with  a  standard 
deviation  of  +4%  around  the  average  displacement,  2^^^.  Test 
results  also  verify  the  simple  constant-C^  model  from  eq. 
2-4,  eq.  2-5  and  eq.  2-20  where  A  simple  curve-fit 

model  was  proposed  for  use  in  determining  velocity. 
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displacement  and  acceleration  during  startup  in  a  drop  test. 
Experimentally  determined  values  of  hydrodynamic  mass  have 
been  presented  and  compared  to  simple  inviscid  theory,  eq. 
2-4c.  Experimental  results  and  the  constant-C^  model  will  be 
compared  to  more  accurate  digital  computer  model  results  in 
Chap.  4. 
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CHAPTER  4 

DIFFERENTIAL  FORMULATION 

The  purpose  of  this  chapter  is  to  present  the  third  of 
three  studies,  the  numerical  or  differential  formulation  of  a 
cylindrical  body  falling  through  fluid  in  a  closed-end  tube. 
The  results  of  this  numerical  effort  are  compared  to  the 
previous  results  predicted  by  the  control  volume  formulation, 
i.e.,  the  constant— model  from  Chap.  2,  and  the  experimental 
results  of  Chap.  3,  as  well  as  results  obtained  by  other 
researchers . 

Sect.  4.1  Inviscid  Analysis 

To  the  best  of  the  author's  knowledge  there  is  no  known 
closed  form  solution  to  the  exact  mathematical  formulation 
for  a  body  falling  through  an  inviscid  fluid  in  a  closed-end 
tube.  It  was  necessary  then  to  make  a  finite  difference 
inviscid  analysis  to  determine  the  pressure  field  around  the 
body.  Appendix  D  contains  the  details  of  the  inviscid 
analysis.  It  was  anticipated  that  at  the  nose  of  the  body  the 
pressure  field  caused  by  the  frictionless  fluid  would  be  very 
similar  to  the  pressure  field  in  water.  After  carrying  out 
both  investigations  this  result  has  been  verified. 

The  next  sections  discuss  the  viscous  flow  field. 

Sect.  4.2  Factors  Affecting  Viscous  Differential  Simulations 

At  the  most  basic  level,  one  is  very  interested  in  the 
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spatial  resolution  of  both  the  velocity  and  pressure  fields 
surrounding  the  model  as  it  travels  down  the  centerline  of 
the  tube.  Interest  in  computational  methods  derives  from 
their  ability  to  produce  large  quantities  of  detailed 
information  about  laminar  and  turbulent  flows  that  is  not 
readily  available  using  experimental  techniques.  In  addition, 
the  computational  results  can  be  globally  verified 
experimentally  since  total  drag  equals  net  weight  at  terminal 
velocity. 

The  exact  mathematical  formulations  for  conservation  of 
mass,  momentum  and  energy  have  been  known  for  over  100  years. 
The  ultimate  may  be  the,  as  yet  unknown,  exact  solution  for 
these  highly  non-linear  partial  differential  equations.  Far 
short  of  that  we  obtain  approximate  solutions  on  a 
discretized  domain,  a  procedure  that  has  many  theoretical 
and  practical  limitations.  The  present  state  of  the  art  in 
computational  fluid  mechanics  allows  a  reasonably  good  job 
modeling  certain  types  of  turbulent  flows.  A  number  of 
theoretical  and  practical  considerations  are  briefly 
discussed  next. 

Sect.  4.3  Numerical  Approach 

Several  types  of  numerical  simulations  presently  exist 
[41],  Large  Eddy  Simulations  (LES)  model  the  small  eddies  and 
compute  or  resolve  the  large  eddies.  Also  Direct  Numerical 
Simulations  (DNS)  resolve  eddies  of  all  sizes.  The  DNS  models 
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are  limited  to  low  Reynolds  numbers  because  of  limitations  on 
computer  size.  Ref  [36]  points  out  that  DNS  requires  the 
number  of  grid  points  to  be  proportional  to  Reynolds  number 
to  the  2.25  power.  It  is  estimated  in  [36]  that  using  DNS  the 
highest  Reynolds  number  that  can  be  computed  on  a  Cray  2 
machine  is  between  5E3  and  1E4  on  a  256  x  256  x  256  grid.  It 
would  take  about  100  hours  of  Cray  CPU  time  to  get  adequate 
convergence . 

By  comparison,  a  commercially  available  program  called 
FLUENT  [38],  running  on  a  Cray  X-MP/28  with  an  80  x  70  grid, 
to  model  flow  around  the  875  body,  converged  in  2  hours  of 
CPU  at  Re^»5E5.  The  same  case  required  60  CPU  hours  to 
converge  on  a  VAX  8600. 

Neither  the  LES  or  DNS  approach  was  utilized  in  this 
study  primarily  because  of  the  relatively  low  Reynolds 
numbers  that  could  be  examined  in  addition  to  their 
inherently  large  CPU  requirements.  The  approach  used  in  this 
study  was  to  use  FLUENT  to  calculate  the  flow  field  followed 
by  some  post-processing  to  obtain  shear  and  profile  drag  on 
the  body.  A  brief  development  of  the  analytical  background 
behind  FLUENT  is  presented  next. 

Sect.  4.4  Basic  Equations 

It  was  hypothesized  that  the  field  was  steady,  incompressible 
turbulent  flow.  The  continuity  and  momentum  equations  [13]  are 


; 


1  ap 


p  ax^  axj 


au^  au^ 

axj  ax^ 


wh0re  is  the  eddy  viscosity^  the  molecular 
viscosity,  and  is  tensor  notation  for  time  mean  velocity 
in  the  x^  direction.  The  first  term  on  the  left  of  the 
momentum  equation  is  the  convection  term;  the  next  term 
represents  pressure  gradient,  followed  by  the  diffusion  term. 
Body  forces  such  as  gravity  have  been  neglected.  The 
continuity  and  momentum  equations  describe  the  mean  or 
time— ave raged  flow  when  modeled  numerically  over  a  finite 
number  of  control  volumes  used  to  discretize  the  flow  domain. 

Because  the  flow  is  turbulent  it  is  customary  to  decompose 
the  velocity  into  separate  terms  for  the  mean  and  fluctuating 
components.  This  velocity  splitting  technique  is  called 
Reynolds  decomposition.  The  approach  is  to  use  the  averaged 
form  of  the  Navier-Stokes  equations.  These  two  equations  are 
given  in  the  radial,  r  direction  (V^=u),  and  the  axial,  z 

direction  (V  -w),  as  follows: 

z 


9u  8u 

1 

3p 

'3^u  1 

3u 

u  3^u 

—  +  w —  * 

3r  9z 

P 

3r 

Ur^  r 

3r 

2  ^  2 
r^  3z 
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3w  3w 

1 

8p 

'3^w  1 

- -  ^  __ 

3w  3^w 

II 

P 

-  -f.  y) 

3z 

Lr^  r 

■  ^ 

3r  3z^ 

'  8  -  8  — 2  ^  -  ' 

-  —  u'w'  +  —  w'  +  —  u'w' 

^  3r  3z  r  J  4-2 

The  Navier-Stokes  equations  are  similar  to  the  momentum 

equations  except  for  the  addition  of  the  Reynolds  stress 

terms,  given  as  the  last  term  on  the  right-hand  side. 

The  continuity  equation  for  the  axisymmetric  case  is 

3  3 

—  (ru)  +  —  (rw)  =  0  . 

3r  3z 

It  was  assumed  in  this  study  that  the  flow  was  isothermal. 

Sect.  4.5  Turbulence  Model 

Since  the  Reynolds  number  is  large,  a  turbulence  model 
must  be  simultaneously  solved  over  the  discretized  domain. 
The  popular  two-equation  k-c  model  of  turbulence  was  used  to 
achieve  closure  [13].  The  k-e  turbulence  model  is  valid  when 
all  Reynolds  stresses  are  of  the  same  order  (isotropic  eddy 
viscosity).  It  was  initially  hypothesized  that  this  problem 
did  not  have  any  highly  swirling  flows,  so  was  isotropic. 
Excellent  agreement  between  the  numerical  model  and 
experiment  subsequently  verified  this  assumption. 

In  the  k-e  model  the  mean  velocity  is  the  only  velocity 
used.  The  effects  of  turbulence,  contained  in  the  Reynolds 
stress  terms,  are  addressed  through  an  "effective  viscosity" 
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term,  that  is  added  to  the  laminar  or  molecular 
viscosity,  ^1*  *rhe  decomposition  of  viscosity  into  these  two 
components  is  known  as  the  "effective  viscosity  hypothesis.” 
The  implementation  of  this  technique  to  any  region  dominated 
by  inertial  effects  (high  local  Reynolds  number)  leads  to  the 
construction  of  the  turbulent  kinetic  energy  (TKE),  which  is 
represented  by  k,  and  the  dissipation  rate  of  TKE,  called  e. 
The  governing  equations  for  k  and  t  are  models  of  the 
Navier-Stokes  equations  and  are  solved  outside  of  the  viscous 
region.  Rodi  [13]  gives  the  k-e  model  in  tensor  form: 


3(Ujk) 
3x  . 


P 


3k 


/ 


3(U^C) 

t 

__  r’  _ 

D  — . 

a 

ax. 

C 

1 

r  — 

^  k 

ax j  . 

where  P 

/ 

=  ^t 

V 

au. 

ax. 

au.  ' 
+  — ^ 
ax^  , 

3U. 

1 

- .  • 

3x  . 

3 

P  is  the  production  of  k.  Also,  is 
the  Prandtl-Kolmogorov  relation: 


4-3 


related  to  k  and  e  by 


■V 


t 


9 


where  the  empirical  constants  in  the  above  are  given  by 


-  1.44,  C2 


1.92,  =  1.0,  Cg  =  1.3, and  =  0.09. 
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Very  close  to  the  wall  the  flow  is  dominated  by  viscous 
effects.  It  was  assumed  that  in  the  viscous  region  the 
logarithmic  law— of— the— wall  applied.  The  k— e  model  was  used 
only  far  from  the  wall,  where  >>  .  Very  near  the  wall 

the  approach  was  to  use  the  logarithmic  law— of— the— wall  from 
the  wall,  y-0,  out  to  y  «  just  outside  the  viscous 

sub-layer.  The  values  of  k  and  c  obtained  at  y^  are  the 
boundary  conditions  for  eg.  4—3.  At  y^  it  was  assumed  that 
the  production  of  TKE  was  exactly  equal  to  the  dissipation  of 
TKE;  in  other  words  there  was  local  equilibrium.  If  k^  and 
denote  the  TKE  and  dissipation  rate  at  point  c,  then  Rodi  has 
shown  that 


w* 


1/2 


and 


w*' 


Ky, 


1  /2 

where  the  friction  velocity,  w*  •=  (t^p)  '  and  K  »  0.41  is 
von  Karman's  constant.  The  law-of-the-wall  used  in  this  study: 


-  2.5  InOr"^) 

for  w'*’  »  normalized  axial  velocity  «  w/w*  , 
and  r'*’  ■  pyw*/^j^ 

for  y  ■  radial  distance  from  the  wall  out  to  a  near-wall 
point  P. 


Sect.  4.6  Method  of  Solution 

The  differential  equations  are  integrated  over  the 
computational  cells  (the  finite  control  volumes)  that 
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comprise  the  flow  domain.  That  is  why  this  method  is  known  as 
the  "finite  volume"  method.  The  core  of  the  computer  method 
is  the  program  FLUENT  [38],  which  is  based  upon  the  SIMPLE 
algorithm  discussed  by  Patankar  [37].  FLUENT  allows  the 
solution  of  elliptic  partial  differential  equations  of  the 
form 


DIV 

(pV<P)  = 

DIV 

(r^  GRAD*)  +  S^  , 

4-4 

where  the 

first, 

second  and  third  terms  are 

the  convection. 

diffusion 

and  source 

terms  respectively.  In 

this  problem  *, 

the  flow  parameter  of 

interest,  is  u,w,k  or 

c,  while  is 

called  the 

dif fusivity 

constant  and  S^  is  the 

source  term.  It 

is  possible  to  rewrite  eq.  4-2  and  4-3  in 

the  form  of  eq. 

4-4.  Table 

7  summarizes  the  results. 

Table 

7. 

Diffusion  and  Source  Terms 

Equation 

4> 

^<t> 

Continuity 

1 

0 

0 

Navier- 
Stokes  r; 

u 

"t 

9  f  3u']  ^  ®  f  ^ 

—  //. —  + - - 

3z  1  3zJ  r  3r  1  3rj  r  3r 

Navier- 
Stokes  z: 

w 

^t 

3  r  3w^  1  3  f 

—  + - 

3z  L  3zj  r  3r  1,  3rJ 

3p 

3z 

TKE 
k : 

k 

P  -  pe 

Dissipation 

Rate  c:  e 

e 

-  (C^P  -  C2PC) 

K 

Note;  P  and  all  the  constants  have  been  defined  in  Sect.  4.5. 
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Discretization  results  in  a  system  of  algebraic  equations 
written  in  terms  of  the  unknown  properties  (fip,  at  a  point  p  : 

♦n  ^  (^i  -  ^  ^i  '*’i  ■  ' 

Pi  ^  i 

where  the  summation  is  over  the  neighboring  computational 
cells,  i  »  N,S,E,W,F,B.  (i.e.,  North,  South,  East,  West,  Front 
and  Back).  The  coefficients,  ,  contain  terms  from  the 
convective  and  diffusive  fluxes.  The  source  terms,  and  Sp 
are  related  to  through 

+  s  *p- 

In  order  to  interpolate  between  grid  points,  which  are  at 
the  center  of  each  cell,  and  to  calculate  the  derivatives  of 
the  flow  variables,  a  power  law  differencing  scheme  was 
utilized.  The  dependent  variables  (p,k,e,v^)  were  calculated 
and  stored  at  the  center  of  the  cells.  Velocity  components 
were  calculated  and  stored  on  the  cell  boundaries,  an 
arrangement  known  as  a  "staggered  grid." 

Boundary  conditions  for  velocity  require  that  the  velocity 
be  specified  at  all  points  on  the  boundary.  Along  the 
centerline  the  radial  component  of  velocity  is  set  to  zero. 
At  the  other  boundaries:  1)  no-slip  velocity  at  the  tube  wall 
and  on  the  model;  2)  axial  velocity  profile  known  or  assumed 
at  the  inlet  together  with  the  zero-gradient  condition;  3)  at 
the  outlet,  the  zero-gradient  condition  is  applied. 

Solution  of  the  simultaneous  set  of  algebraic  equations 
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was  accomplished  using  a  semi-implicit  iterative  scheme, 
described  in  the  next  section. 


Sect.  4.6.1  The  Pressure-Correction  Scheme 

It  is  convenient  to  derive  an  equation  that  allows 
pressure  to  be  calculated  directly  from  the  velocity  field. 


The  continuity  equation  is  combined  with  the 

divergence 

of 

the  momentum 

equation  and  is  rewritten 

into  a 

Poisson-type 

equation : 

3^p  - 

9U^  9U^ 

4 

-5 

”=  —  p 

9x j  9x j 

9Xj  9x^ 

Once  eq. 

4-5  is  integrated 

ove  r  a 

control  volume 

and 

subsequently 

discretized  it  is 

called 

the  p' 

equation 

or 

pressure  correction  equation. 

It  is  , 

then 

possible 

to 

calculate  p'  at  a  point  P  given  the  velocities  at  neighboring 

nodes.  The  FLUENT  solution  algorithm  is  as  follows: 

* 

1.  Guess  a  pressure  field  p  . 

*  * 

2.  Solve  the  Navier-Stokes  equations  to  obtain  guesses  u  ,w  . 

3.  Solve  the  p'  equation. 

* 

4.  Update  pressure,  p  »  p  +  P' 

5.  Calculate  velocity  corrections,  u'  and  w'  from  p' 

*  * 

6.  Update  velocity;  u»=u  +u',w  =  w  +w' 

7.  Solve  the  k  and  t  equations  with  the  updated  velocities. 

8.  Directly  calculate  the  field  from  the  k  and  e  fields. 

9.  Let  p  become  the  new  guess  for  pressure  and  start  over. 
Steps  2-9  are  repeated  until  the  pressure  correction  term 
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has  decreased  to  a  required  value  for  each  point  on  the  grid. 
At  this  point  local  continuity  has  been  satisfied.  Global 
continuity  is  checked  by  integrating  the  velocity  over  the 
outlet  boundary  condition. 

Sect.  4.7  Discussion  of  Computer  Simulation 

The  author's  experience  with  the  computer  algorithm 
indicates  that  the  velocity  field  settles  down  to  essentially 
the  converged  value  after  only  200—400  iterations,  regardless 
of  Reynolds  number.  The  pressure  field  converges  much  slower, 
taking  up  to  8000  iterations  in  a  number  of  cases,  the  number 
of  iterations  increasing  with  increasing  Reynolds  number. 

Since  the  convergence  and  accuracy  of  the  numerical 
simulation  depend  on  grid  size  and  the  location  of  the  inlet 
and  outlet  planes,  these  effects  were  studied  during  this 


investigation . 

Several 

grid  layouts  were 

used  during 

the 

course 

of  the 

study. 

The  location 

of 

the 

inlet  plane 

was 

varied 

with 

very  little  effect 

on 

the 

flow  field 

or 

calculated  drag  force.  It  was  also  determined,  by  numerical 
experimentation,  that  drag  force  was  relatively  insensitive 
to  the  location  of  the  outlet  plane.  Grid  size  was  also 
varied.  Table  8  presents  results  from  the  grid  size  and 
inlet  and  outlet  locations  and  the  effect  on  drag  for  the  875 
model.  The  875  body  was  selected  since  it  had  the  smallest 
annular  9ciPr  causing  the  flow  to  distort  a  maximum  amount  to 
go  from  the  tube  into  the  annulus,  in  each  case  presented  in 
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Table 


it  was  assumed  that  the  875  model  was  traveling  at 


5.1E5. 


a  velocity  of  75  ft/s  in  water,  corresponding  to  Re^ 

Overall  pressure  drag  is  represented  by  and  total  shear 
drag  by  D^,  both  given  in  pounds. 

Table  8.  Grid  Size  and  Inlet/Outlet  Location  Effects  on  Drag 


Case 

Grid  Size  Inlet  Location  Outlet  Location 

Axial  X  from  the  nose  from  the  tail 

Radial  in.  m.  p,  lb. 

lb. 

1 

80X70  1.0 

12.5  50.88 

2.173 

2 

15. o' 

50.99 

2.172 

3 

80X70  12.5 

1.0  50.87 

2.173 

4 

5.0  50.97 

2.175 

5 

15.0  50.99 

2.172 

6 

100X100  15.0 

15.0  50.95 

2.169 

7 

80X70 

50.99 

2.172 

8 

70X50 

51.29 

2.194 

8 

58X50 

51.37 

2.200 

9 

40X35 

53.54 

2.362 

10 

30X26 

63.53 

2.509 

Comparing  cases  1  and  2 

it  is  apparent  that  the  inlet 

plane 

can  be  located  anywhere 

from  1"  to  15"  from 

the  nose 

with  : 

less  than  a  0.25%  change 

in  drag.  Location  of  the  outlet 

plane 

was  studied  under  cases 

3,  4  and  5,  with  the  similar 

result  that  drag  is  insensitive  to  outlet  plane  location.  The 
author  suggests  that  the  findings  presented  in  cases  1 
through  5  are  an  important  result  of  this  study.  The 
conclusion  that  drag  is  independent  of  the  location  of  inlet 
and  outlet  planes  must  be  qualified  by  observing  that  this 
study  was  concerned  only  with  highly  streamlined  bodies  and 
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that  a  bluff  body  could  not  be  expected  to  behave  in  a 
similar  fashion. 

Grid  size  was  studied  in  cases  6  through  10.  The  inlet 
plane  was  located  15"  ahead  of  the  nose  of  the  model.  The 
outlet  plane  was  placed  15"  downstream  from  the  extreme  end 
of  the  tail  to  capture  the  wake.  Originally,  in  order  to 

establish  a  baseline,  the  finest  grid  allowed  by  FLUENT, 

figure  22  ,  having  100  divisions  in  the  radial  and  100 
divisions  in  the  axial  direction,  was  run  on  the  875  model. 
The  coarse  30X26  grid  in  case  10,  over— predicts  drag  by  about 
24%  compared  to  the  finest  100X100  grid  in  case  6.  Since  case 
8  deviated  by  only  about  0.5%  from  case  6,  with  65%  fewer 
grid  points  and  requiring  less  than  10%  of  the  CPU  time,  it 

was  decided  to  use  the  70X50  grid  for  all  the  Reynolds 

numbers  that  were  studied.  The  70X50  grid  required  between  24 
and  48  CPU  hours  to  converge  on  the  VAX  8600,  depending  on 
the  Reynolds  number.  The  larger  the  Reynolds  number  the 
slower  the  pressure  field  converges. 

Sect.  4.8  Presentation  and  Discussion  of  Results 

The  625  model  was  studied  in  the  most  detail  primarily 

3 

because  the  dimensionless  groups  n2*A/D  «1.64  and 
Hg-p-D/d-l . 6 ,  Sect.  3.7,  were  of  the  most  immediate  interest 
to  the  Navy.  This  case  represented  a  prototype  21"  diameter 
torpedo  inside  a  27"  diameter  torpedo  tube,  a  realistic 
combination. 
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As  discussed  in  Sect.  3.10,  the  625  body  also  had  the  most 
interesting  acceleration  region,  being  the  greatest  in 
magnitude  and  over  the  longest  time,  of  the  three  models.  The 
625  also  had  the  highest  terminal  velocity,  56  in./s,  and  the 
highest  gap  Reynolds  number,  ^6^=  22200,  using  figure  12  : 

2 

wD  P 

Re  “  — —  !  where  w  =  —x —  V  and  D  =  D-d; 

5  V  e^-i 

s' 

then.  Re  *=  -  Re,  . 

9  S+1 

For  these  reasons  it  was  decided  that  the  625  would 
receive  the  most  detailed  study  of  the  three  configurations. 
The  750  and  875  were  also  studied  for  lElb^<  Re^  <5E10^  in 
order  to  determine  the  functional  dependence  of  in  terms 

of  p  and  Re  . 

9 

There  was  an  assumption  made  in  Sect.  2.10  that  nose  drag 
and  tail  drag  are  independent  of  the  length  of  the  annular 
region.  The  exact  hypothesis  in  Sect.  2.10  was  that  total 
drag  can  be  found  by  independently  calculating  nose  drag, 
annular  drag  and  tail  drag  and  finally  adding  them  together. 
Of  course  the  question  is  whether  the  actual  composite  model 
behaves  in  this  ideal  independent  fashion.  To  test  this 
hypothesis  a  special  case  was  run  on  FLUENT,  called  625 
Shorty.  Shorty,  figure  23  ,  was  basically  a  nose  and  tail, 
with  no  annular  region.  Comparison  of  pressure  coefficients 
over  a  wide  range  of  with  the  full  625  model  revealed 
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that  the  hypothesis  was  valid. 

A  graph  of  Cp  vs.  axial  position  for  various  values  of 
turbulent  Shorty  is  presented  in  figure  24  .  The 

two  regions  of  interest,  the  nose  and  tail,  are  presented  and 
discussed  in  the  following. 

Figure  25  shows  a  composite  overlay  of  Cp  on  the  nose. 
The  inviscid  prediction  computed  by  finite  differences  is 
compared  to  FLUENT  predictions  for  a  wide  range  of  turbulent 
Reynolds  numbers.  The  overall  prediction  of  nose  drag 
compares  very  favorably,  showing  independence  from  Reynolds 
number  and  only  about  +4%  maximum  variation  between  the 
inviscid  and  FLUENT  at  the  back  of  the  nose.  In  addition  to 
corresponding  well  with  the  full  625  model,  the  nose  drag  on 
Shorty  compares  very  well  with  the  simple  formula,  eq.  2-9, 
where  it  was  assumed  that  there  was  essentially  no  pressure 
drop  over  the  nose  due  to  friction.  This  important  result 
conveniently  allows  use  of  eq.  2-9  with  K^*Nose  loss 

coefficient-0,  as  in  eq.  2-20,  to  compute  the  nose  drag 

coefficient, 

Figure  26  is  a  composite  of  various  turbulent  Reynolds 
numbers  for  Cp  vs.  axial  position  on  the  tail  of  Shorty.  The 
tail  curves  from  figure  24  have  simply  been  shifted  up  or 

down  by  a  constant  in  order  to  show  that  C^  from  eq.  2-11, 

with  K^-Tail  loss  coefficient-1.0,  as  in  eq.  2-20  accurately 
predicts  the  tail  drag.  This  is  another  important  result,  and 
verifies  that  tail  drag  coefficient,  C^,  may  be  found  using 
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the  constant-C^  model  of  eq.  2-20. 

Figure  27  presents  Cp  vs.  axial  position  for  the  full 
625  model  at  various  values  of  comparison  the 
inviscid  prediction  was  plotted  with  the  FLUENT  results.  The 
nose  is  seen  to  behave  inviscidly,  as  in  Shorty.  The  straight 
section  then  opens  up  into  a  series  of  straight  lines  of 
differing  slopes,  and  the  tail  curves  look  similar  except 
that  they  are  shifted  by  a  constant.  Each  of  the  three 
regions  of  interest,  nose,  straight  and  tail  are  discussed 
next. 

Figure  28  contains  a  composite  overlay  of  Cp  on  the  nose 
of  the  full  625  model  for  the  inviscid  and  various  turbulent 
Reynolds  numbers.  The  inviscid  prediction  compares  very 
favorably  to  the  FLUENT  results,  as  in  figure  25  .  There  is 
especially  excellent  correlation  at  the  back  of  the  nose, 
resulting  in  about  +1%  maximum  deviation.  This  presents  the 
final  piece  of  evidence  that  eq.  2—20  can  be  used  to  compute 
the  nose  drag  coefficient, 

Figure  29  overlays  Cp  for  the  tail  of  the  625  body, 
showing  that  the  tails  are  in  fact  shifted  by  a  constant,  as 
in  the  Shorty  tail  results  of  figure  26  .  Comparing  the  tail 
results  from  Shorty  in  figure  26  to  the  tail  results  of 
the  full  625  model  in  figure  29  shows  that  tail  drag  varies 
only  about  ±2%  between  the  two.  This  is  an  extremely 
important  result  since  it  verifies  the  hypothesis  that 
tail  drag  is  essentially  unaffected  by  the  length  of  the 
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straight  section. 

NOW  that  the  nose  and  tail  sections  have  been  discussed 
and  it  has  been  shown  that  they  both  behave  as  hypothesized 
in  Chap.  2,  it  is  time  to  discuss  straight  section  drag. 
Verification  of  Kotlow's  [10]  straight  section  model  is 
obtained  from  figure  27  .  Using  Kotlow's  formulas  to  compute 
C  -Straight  drag  coefficient,  eq.  2-17,  the  FLUENT  results 
are  predicted  by  +5%,  which  is  considered  very  good.  Kotlow 
points  out  that  neglect  of  the  developing  entrance  region 
will  produce  about  +4%  error  in  drag.  Observe  in  figure  27 
that  there  is  essentially  no  discernible  change  in  slope  in 
the  drag  coefficient  computed  by  FLUENT  along  the  straight 
section.  Kotlow's  assertion  that,  for  purposes  of  computation 
of  drag,  the  entire  straight  region  may  be  treated  as  fully 
developed  is  then  verified  by  this  study. 

While  it  is  apparent  that  the  velocity  profiles  will  be 
very  different  along  the  straight  section,  it  is  not  apparent 
that  the  drag  coefficient  is  essentially  independent  of  axial 
position  in  the  straight  section.  The  next  paragraphs  discuss 
some  FLUENT  results  for  velocity  profiles  at  various 
locations,  called  'stations',  in  the  flow  field. 

Recall  that  in  a  coordinate  system  fixed  to  the  model,  the 
625  body  has  a  zero  velocity  and  the  tube  wall/inlet  plane 
have  a  velocity  of  4.67  ft/s,  which  corresponds  to  22200 

at  terminal  velocity  in  a  1"  tube  filled  with  water.  Figure 
30  was  run  on  FLUENT  for  the  terminal  velocity  case  on  the 
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100X100  grid  and  shows  velocity  profiles  just  ahead  of  the 
nose  (station  13),  at  the  beginning  of  the  straight  annulus 
(station  25)  and  at  the  middle  of  the  straight  annulus 
(station  40).  The  flow  is  almost  a  slug  flow  as  it  enters  the 
straight  section  and  rapidly  begins  to  form  into  a  classic 
curved  profile  due  to  friction  with  the  wall  and  model. 

Figure  31  is  an  extension  of  figure  30  displaying 
velocity  profiles  at  the  end  of  the  annulus  (station  56), 
middle  of  the  tail  (station  62),  just  behind  the  model 
(station  74)  and  at  the  outlet  plane  positioned  15  inches 
downstream  of  the  tail  (station  100).  The  flow  is  seen  to 
reach  a  nicely  developed  profile  in  the  straight  section  and 
then  expand  out  into  the  tail  region.  Baclt  out  in  the  tube, 
at  station  100,  the  effects  of  the  high  gap  velocity  are 
fading  but  still  noticeable. 

A  comparison  of  velocity  profiles  with  the  predictions  of 
Kotlow  [10]  is  also  of  interest.  Figure  32  shows  Kotlow's 
fully  developed  profile  compared  to  the  profile  predicted  by 
the  differential  model.  The  maximum  difference  is  about  +3.5% 
between  the  two  velocity  profiles,  which  is  considered  quite 
good. 

Figure  33  presents  the  ratio  of  profile  drag  to  total 
drag  for  the  three  bodies.  For  the  625  body,  the  profile  drag 
is  seen  to  start  at  80%  for  ROg^lO  and  increase  to  85%  for 
Reg»5ElO®.  Figure  33  highlights  the  fact  that  drag  is 
primarily  caused  by  pressure  forces. 
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Drag  coefficients 

figure  13  .  A  curve 
and  lE5<Reg<5E6,  was 

are  plotted  vs.  gap  Reynolds  number 

fit  in  the  range  of  interest,  1.14<8<1 

developed  to  predict  drag  coefficient: 

Define; 

\ 

A 

and 

Re  =  -  Re  .  • 

^  2 

A-a 

9  B.1 

Then: 

exp(55//Reg 

)  Reg'^, 

4-7 

where 

b  -  0.01 

O 

• 

O 

+ 

_l)-0-6  , 

and 

r  -  2.92 

-  0.8536 

(yields  +5% 

accuracy  for  C^) 

or 

r  -  -0.60698^  +  0. 

76786  +  1.8652  (+3%  for  C^). 

625 

750 

875 

1.600 

1. 

333 

1.143 

X 

1.641 

2. 

286 

4.267 

b 

0.02828 

0. 

02204 

0.01688 

r 

1.54 

1. 

81 

1.95 
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It  is  of  interest  to  note  in  figure  13  that  if  a  power 
law  relation  is  assumed, 

C,  =  k  Re  (k  and  n  are  constants), 

d  g 

the  best  fit  values  in  the  range  10^<Reg<5El0®  are 


625 

750 

875 

& 

1.600 

1.333 

1.143 

k 

3.833 

31.797 

998.15 

n 

0.06475 

0.1458 

0.2791 

It  is  clear  from  this  table  that  is  a  strong  function  of 
3,  as  discussed  next. 

When  3  is  small,  the  gap  clearance  is  small,  with  the  result 
that  becomes  large  and  highly  sensitive  to  changes  in  body 
velocity.  The  reason  is  that  the  gap  velocity  is  greatly 
magnified  by  the  small  clearance  area  compared  to  the  tube 
area.  Hence, for  small  values  of  3/  is  very  sensitive  to 
changes  in  Re^. 

However,  when  3  is  large,  the  gap  area  becomes  large 
compared  to  the  tube  area  and  the  gap  velocity  is  less 
sensitive  to  changes  in  body  velocity.  Hence, for  large  values 
of  3,  is  less  sensitive  to  changes  in  Re^. 
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Summary 


This  chapter  has  presented  and  discussed  drag  computed  by 
finite  differences  to  both  the  inviscid  and  viscous  problem 
of  a  body  falling  through  fluid  in  a  closed-end  tube.  It  has 
been  shown  that  nose  drag  in  an  inviscid  flow  is  the  same  as 
nose  drag  in  viscous  flow,  i.e.  the  effect  of  shear  on  the 
nose  is  negligible.  It  has  also  been  shown  that  nose  drag  and 
tail  drag  are  independent  of  the  length  of  the  cylindrical 
section.  This  conclusion  allows  for  convenient  independent 
calculation  of  nose,  straight  and  tail  section  drag,  eq. 
2-19.  Total  drag  is  then  the  sum  of  the  independently 
computed  components.  It  has  been  further  shown  that  nose  and 
tail  drag  for  turbulent  flow  are  independent  of  Reynolds 
number.  Drag  on  the  straight  section  does  depend  on  Reynolds 
number  in  accordance  with  Kotlow's  prediction,  eq.  2-16  and 
2-17.  Finally  drag  coefficients  have  been  plotted  as  a 
function  of  P“D/d  and  Reynolds  number  using  both  the 
constant— model  of  Chap.  2  and  the  numerical  simulation  of 
Chap.  4.  A  curve-fit  has  been  proposed  for  the  range  of  P  and 
Reynolds  numbers  of  interest  to  the  Navy. 
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CONCLUSIONS  AND  RECOMMENDATIONS 

Conclusions 

1.  This  study  has  shown  that  nose  drag  and  tail  drag  are 
independent  of  the  length  of  the  straight  annular  region. 
This  conclusion  is  specific  to  a  model  having  a  hemispherical 
nose,  parallel  midsection  and  conical  tail,  moving 
concentrically  in  a  fluid-filled  tube.  This  conclusion  was 
confirmed  numerically  for  1.14<(5<1.6  by  comparison  of  drag 
from  calculations  on  a  body,  called  SHORTY,  having 
essentially  zero  length  in  the  cylindical  section,  with  drag 
calculations  on  bodies  of  finite  length. 

2.  Drag  coefficients  predicted  by  control  volume  methods, 
experimental  drop  tests,  and  numerical  simulation  were 
consistent  (within  +5%)  for  turbulent  Reynolds  numbers  in  the 
range  lE5<Reg<5E6.  Estimates  obtained  in  this  study  are 
consistent  (within  +5%)  with  those  of  previous  investigators. 

3.  A  curve-fit  has  been  proposed  to  predict  drag 
coefficients  for  bodies  in  the  range  1.14<e<1.6  with 
lE5<Reg<5E6,  which  has  +5%  accuracy.  This  is  the  range  of 
most  interest  to  the  Navy  for  torpedo  launching. 

4.  The  variable-C^  analysis  is  only  approximately  correct. 
The  physics  of  the  accelerating  flow  around  the  model  are  not 
well  understood.  Differential  analysis  of  the  unsteady 
problem,  where  the  vehicle  accelerates  from  rest  to  the 
terminal  velocity,,  was  not  attempted  in  this  study. 
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Recommendations 

1.  Future  investigators  may  want  to  attempt  a  numerical 

solution  of  the  unsteady  case. 

2.  A  companion  problem  to  the  one  discussed  in  this  study 
involves  the  investigation  of  a  model  in  a  finite  length 
open-end  tube.  In  this  case  the  pressure  drop  developed 
across  the  model  exactly  equals  the  pressure  drop  across  the 
tube.  Flow  is  thereby  induced  in  the  tube  in  the  far  field.  A 
numerical  solution  of  the  flow  field  and  model  drag  may  be 
attempted.  The  one  dimensional  constant-C^  analysis  developed 
in  this  study  could  be  generalized  to  handle  this  case. 

3.  The  problem  described  in  2.  above  could  be  examined 
for  the  unsteady  case  also. 

4.  Another  companion  problem,  of  interest  to  the  Navy,  is 
the  case  of  a  favorable  pressure  gradient ,  i.e. ,  from  a  pump, 
applied  to  the  model  in  order  to  push  it  out  of  the  tube. 
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Figure  1.  Basic  Torpedo  Launch  Configurations 
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Figure  6b.  Converging-Diverging  Nozzle,  of  area  Ag,  used  to 
model  the  tube-model  flow  in  figure  (6a) 


Figure  6.  Converging-Diverging  Nozzle  model  of  actual 
Tube-Model  flow  (used  to  model  compressible 
effects,  as  discussed  in  Appendix  A) 
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Figure  9.  Dist.,  Z,  vs.  Time,  t,  for  the  625  model,  Averaged  Experimental  Data 
compared  to  Constant-C,  model,  Variable-C,  model  and  Curve-Fit  model 
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Figure  11.  Distance,  Z,vs.  Time,  t,  for  the  875  model 
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Figure  15.  Displacement,  Velocity  and  Acceleration  vs.  Time  for  the  625  model 
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Figure  16.  Displacement,  Velocity  and  Acceleration  vs. 
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Figure  17.  Displacement,  Velocity  and  Acceleration  vs.  Time  for  the  875  model 
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Figure  18.  Drag  Coefficient  vs.  Time,  predicted 


Az  =  0.03"  Uniform  Grid 

Ar  =  0.005"  625  model 

Numbers  inside  the  g-rid  signify  node  numbers;  for  example, 
node  62  in  the  r  direction  is  the  outer  radius  of  the  model, 
corresponding  to  a  physical  distance  of  62*0.005  =  0.031". 


Figure  19.  Inviscid  Grid  Layout  for  the  625  model 
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Figure  21  Inviscid  Pressure  Coefficient  vs.  Axial  Position  for  the  625  model 
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Figure  22.  FLUENT  grid  layout  for  the  625  model 


Numbers  inside  the  grid  signify  node  numbers 
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Figure  24.  Pressure 
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station  numbers  correlate  to  axial  nodes,  see  figure  22 
All  curves  computed  using  FLUENT,  Section  4.8 
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Figure  30.  Velocity  Profiles,  Stations  13,  25  &  40,  for  the  625  model 
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Figure  33.  Profile  Drag/Total  Drag  Ratio  for  the  625,  750  &  875 


APPENDIX  A 


FLUID  PROPERTIES  AND  EXPERIMENTAL  CONSIDERATIONS 

This  appendix  presents  several  tables  of  viscosities  for 
various  fluids  and  discusses  experimental  methods  of 
achieving  low  viscosities.  Several  other  experimental  issues 
are  briefly  discussed  including  compressibility  effects  and 
cavitation. 

Gaseous  Medium 

In  order  to  maximize  the  Reynolds  number,  a  fluid  with 
very  low  kinematic  viscosity  is  desirable.  In  this 

connection,  it  is  instructive  to  examine  a  relation  from 
kinetic  theory.  White  [40]  for  kinematic  viscosity  of  gas 
(there  is  no  comparable  formula  for  liquids): 


1.442E-5 

t3/2 

ft^ 

2 

a  2  p 

sec 

(T  V 

-0.145 

(T 

B  =  1.147 


—  +0.5 


-2.0 


where  a  »  collision  diameter.  Angstroms;  M  =  molecular  wt; 

T  «  absolute  temperature,  “R;  T  »=  effective  temperature,  K; 
p  *  pressure,  psi;  2  *  collision  integral,  dimensionless. 

At  standard  conditions,  p  *  14.7  psi,  T  =  530  ®R,  we  have: 


131 


(T  =  294  K  in  2  equation) 


Table  A-1.  Kinematic  Viscosities  from  Kinetic  Theory 


Gas 

M 

a 

52 

Ai  r 

29 

3.711 

78.6 

1.003 

1.61E-4 

Ar 

39.9 

3.542 

93.3 

1.046 

1.44E-4 

CI2 

70.9 

4.217 

316.0 

1.648 

4.85E-5 

CH4 

16 

3.758 

148.6 

1.202 

1.76E-4 

CO2 

44 

3.941 

195.2 

1.329 

8.74E-5 

CO 

28 

3.690 

91.7 

1.042 

1.60E-4 

«2 

2 

2.827 

59.7 

0.944 

1.13E-3 

He 

4 

2.551 

10.22 

0.706 

1.30E-3 

^2 

28 

3.798 

71.4 

0.981 

2.26E-4 

NO 

30 

3.492 

116.7 

1.113 

1.61E-4 

N2O 

44 

3.828 

232.4 

1.430 

8.61E-5 

O2 

32 

3.467 

106.7 

1.085 

1.62E-4 

Thus 

from  the  above  table  it  is 

seen  that  chlorine 

offers 

the 

lowest  viscosity.  Indeed  it  is  difficult  to  get 

below 

V-2E 

-5 

for  any 

gas  at  standard  pressure  and 

temperature , 

given  the  range  of  a 

and  T  .  Chlorine  under  1000 
e 

atm  pressure  at  room  temperature  (294  K)  gives  v=5E-8.  On  the 
other  hand  it  could  be  done  with  chlorine  at  1  atm.  pressure 
and  liquid  helium  temperatures  (3  K).  From  the  formula,  one 
should  look  for  a  gas  with  large  collision  diameter,  c,  large 
collision  integral,  2,  and  a  large  molecular  weight,  M. 

In  the  event  that  a  gas  is  chosen  as  the  fluid  medium, 
one  must  attempt  to  match  both  the  Reynolds  and  Mach  numbers. 
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The  point  in  the  flow  field  that  will  experience  the  highest 
velocity  (and  therefore  greatest  Mach  number)  occurs  at  that 
cross-section  of  the  body  having  the  greatest  diameter.  This 
study  was  concerned  with  just  the  incompressible  case,  so 
velocity  must  be  kept  well  below  the  speed  of  sound  of  the 
fluid.  An  approximate  analysis  can  be  performed  to  relate 
Reynolds  number  as  a  function  of  3  and  the  physical  constants 
of  the  gas.  The  tube-model  interaction  is  similar,  figure  6 
to  a  one  dimensional  analysis  of  a  converging-diverging 
nozzle  with  Mach  number  *  0.2  at  the  throat.  To  push  the 
allowable  Reynolds  number  even  higher  demands  that  a  gas  be 
found  having  a  combination  of  the  highest  speed  of  sound,  c, 
and  lowest  kinematic  viscosity. 


V  =  0.2c 

g 


3' 


-  1 


vd  D 

V,  Re  ,  =  — ,  d  -  R 


49700 


M 


3.716  D  (3^-1) 


Re 


d,max 


IM  J 


where  D*Tube  diameter,  inches;  M^molecular  weight;  T»absolute 
temperature,  "R;  v-kinematic  viscosity,  ft  /s.  Specifically, 
a  gas  should  have  a  low  molecular  weight,  low  kinematic 
viscosity  and  a  high  specific  heat  ratio.  A  careful  search 
was  made  [42],  [43],  [44]  and  the  results  for  37  different 
gases  follow:  (with  T  -  530‘'R,  D  =  1  in.,  3  “  1.6) 
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Table 


A-2.  Kinematic  Viscosities  and  Limiting  Reynolds 


Numbers  for  Various  Gases 


Gas 

Formula 

M 

Y 

V  *10^ 

Re 

*10^ 

Acetylene 

^2«2 

26 

1.24 

100. 

71 

Ai  r 

29 

1.40 

163. 

44 

Ammonia 

NH3 

17 

1.29 

155. 

58 

Argon 

Ar 

39.9 

1.66 

131. 

51 

n-Butane 

^4«10 

58.1 

1.10 

31.2 

143 

Carbon  dioxide 

CO2 

44 

1.28 

82.7 

68 

Carbon  monoxide 

CO 

28 

1.40 

169. 

43 

Chlorine 

C12 

70.9 

1.33 

51.9 

87 

Deuterium 

^2 

4 

1.40 

1700. 

16 

Ethane 

^2^6 

30.1 

1.19 

78.9 

82 

Ethylene 

^2»4 

28.1 

1.24 

93.4 

73 

Fluorine 

^2 

38 

1.36 

166. 

37 

R-11^ 

CCI3F 

137.4 

1.14 

170. 

(  1) 

17 

16.5 

180 

R-12  CCI2F2 

120.9 

1.14 

20.7 

(107) 

153 

30.5 

104 

R-22  CHCIF2 

86.5 

1.19 

18.3 

(  22) 

209 

41.5 

92 

R-113  CCl 

jF-CClFj 

187.4 

1.12 

318. 

(  0) 

8 

25.5 

99 

1.  Two  values 

are  given 

for 

each  R-XYZ 

.  The 

first  is 

for 

saturated  vapor  at  -27®F  with  the  approximate  pressure  shown 
in  parentheses.  The  second  value  is  for  superheated 
vapor  at  86“F  and  1  atm  pressure.  All  pressures  (  ),  psia. 
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Table  A-2.  continued 


R-114  C 

2CI2F4 

170.9 

1.09 

78.7 

( 

3). 

33 

16.7 

156 

R-502  CHCIF 

2CCIF2CF3 

111.6 

1.14 

13.2 

( 

31) 

249 

24.5 

134 

R-717 

NH3 

17 

1.29 

104. 

( 

16) 

86 

140. 

64 

Helium 

He 

4 

1.67 

13100. 

16 

Hydrogen 

«2 

2 

1.41 

1170. 

23 

Hydr.  chloride 

HCl 

36.5 

1.39 

108. 

59 

Hydr.  sulfide 

H2S 

34.1 

1.34 

101. 

64 

Krypton 

Kr 

83.8 

1.67 

78.5 

59 

Methane 

CH4 

16 

1.30 

179. 

52 

Methyl  chloride 

CH3CI 

50.5 

1.24 

54.7 

93 

Neon 

Ne 

20.2 

1.67 

417. 

22 

Nitric  oxide 

NO 

30 

1.39 

167. 

42 

Nitrogen 

^2 

28 

1.40 

169. 

43 

Nitrous  oxide 

N2O 

44 

1.27 

89.9 

62 

Nonane 

^3®20 

128.3 

1.15 

20. 

154 

Oxygen 

°2 

32 

1.39 

164. 

42 

Ozone 

O3 

48 

1.27 

71.3 

74 

Propane 

C3«8 

44.1 

1.13 

47.2 

110 

Propylene 

42.1 

1.15 

531. 

10 

Sulfur  dioxide 

SO2 

64.1 

1.39 

53.4 

90 

Xenon 

Xe 

131.3 

1.15 

46.2 

66 

Looking  at 

the  preceding  table 

,  one 

concludes 

that 

it  is 

very  difficult  to  get  a  Reynolds  number  over  200,000  at 
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standard  pressure  and  temperature  without  introducing 
compressible  effects.  This  survey  shows  that  only  marginal 
improvements  can  be  expected  by  varying  the  type  of  gas, 
certainly  not  the  3  orders  of  magnitude  improvement  that  is 
required  to  match  the  prototype  Reynolds  number. 

The  previous  arguments  from  kinetic  theory  can  be  invoked 

to  achieve  a  combination  of  pressure  and  temperature  that 

would  get  the  low  kinematic  viscosity  required.  One 

interesting  alternative  would  be  to  use  room  temperature 

2 

Freon-11  at  about  500  atm  pressure.  This  would  cause  the 
viscosity  to  come  down  to  about  5E-8  and  the  density  would 
simultaneously  increase  on  the  order  of  500  times,  causing 
the  model  to  fall  on  the  order  of  10  ft/sec  at  terminal 
velocity.  Thus  one  would  match  both  the  Reynolds  and  Mach 
numbers  of  the  prototype  precisely,  and  still  have  a 
reasonably  short  tube.  However,  the  practical  difficulties  of 
performing  such  a  test  are  impressive. 

One  is  left  to  conclude  that  at  standard  pressure  and 
temperature,  assuming  a  drop  test  using  any  gas,  there 
results  a  very  long  apparatus  to  achieve  terminal  velocity 
and  that  compressibility  effects  are  unavoidable.  Due  to 
these  concerns,  the  decision  was  made  to  use  a  liquid  in  the 
drop  test  apparatus. 


2.  It  must  be  borne  in  mind  that  the  critical  point  of  F-11 
is  at  388. 4®F,  635  psia,  so  that  at  500  atm.  the  perfect  gas 
law  will  not  apply,  so  the  formula  for  maximum  Reynolds 
number  should  not  be  used  since  it  was  assumed  that  the  speed 
of  sound  was  equal  to  that  of  an  ideal  gas. 
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Liquid  Medium 

Once  again,  the  goal  of  high  Reynolds  numbers  motivated  a 


thorough  search  [42] 

,  [43],  [44] 

for 

liquids 

with  low 

kinematic  viscosities, 

the  results  for 

52 

liquids 

follow: 

(Pressure  =  1  atm)^ 

Table  A-3.  Kinematic  Viscosities 

of 

Various  Liquids 

Liquid 

Formula 

Temperature 

V  *10^ 

“C 

ft*ft/s 

Acetic  Acid 

'=2«4°2 

25 

119 

Acetone 

C3H6O 

40 

43.3 

Alcohol,  butyl 

C^HgOH 

30 

306 

Alcohol,  ethyl 

C2H5OH 

30 

150 

Alcohol,  methyl 

CH2OH 

25 

76.6 

Alcohol,  propyl 

C3H8O 

30 

258 

Ammonia  (~171.) 

NHg 

30 

32.3 

Benzene 

25 

74 

Bismuth 

Bi 

285 

17.3 

Bi 

600 

10.7 

Bromine 

Br2 

20 

34.9 

Carbolic  acid  (phenol) 

25 

803 

Carbon  disulfide 

CS2 

25 

30.7 

Carbon  tetrachloride 

CCI4 

25 

61.8 

Castor  Oil 

9  m  9 

25 

73200 

Chloroform 

CHCI3 

25 

39 

3.  Pressures  other  than  1  atm.  are  shown  in  parenthesis  and 
are  in  psia.  ex.  (~171.)  means  'about  171  psi  absolute'. 
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Table  A- 3*  continued 


Decane 

^10^22 

25 

127 

Dodecane 

^12^26 

25 

196 

Ether 

''4»10° 

25 

33.' 

Ethyl  Acetate 

C^Cl3H502 

25 

33.! 

Ethyl  Bromide 

C2H5Br 

30 

26.: 

Fuel  oil  no.  1^ 

•  •  • 

38 

151 

Fuel  oil  no.  2 

•  •  • 

38 

215 

Fuel  oil  no.  4 

•  •  • 

38 

624 

Fuel  oil  no.  5  (light) 

•  •  • 

38 

3443 

Fuel  oil  nO.  5  (heavy) 

•  «  • 

38 

8070 

Gasoline 

•  •  • 

25 

47 

Glycerin 

'^3»8°3 

25 

81200 

Heptane 

'^7«16 

25 

62.7 

Hexane 

25 

48.8 

Kerosene 

•  •  • 

25 

240 

Lead 

Pb 

350 

26.1 

Pb 

844 

12.0 

Linseed  oil 

•  •  • 

25 

3840 

Mercury 

Hg 

25 

12.1 

Octane 

25 

78.6 

Pentane 

20 

41.3 

Propane 

''3«8 

-45 

20.2 

Propylene 

‘^3‘*6 

-50 

15.9 

Propylene  glycol 

''3''8°2 

25 

4350 

4.  Fuel  oils  are  taken  at  the  minimum  allowed  for  that  grade. 


138 


Table  A- 3.  continued 


R-ll^ 

(  1) 

CCI3F 

1 

0 

42.2 

(  18) 

se^F 

29.6 

R-12 

{  14) 

CCI2F2 

-27“F 

24.1 

(107) 

86“F 

21.1 

R-22 

(  22) 

CHCIF2 

-27“F 

23.1 

(171) 

86°F 

21.1 

R-113 

(  0) 

CCI2F-CCIF2 

-27°F 

81.9 

(  5) 

86“F 

43.7 

R-114 

(  3) 

-27"F 

41.0 

(  37) 

86‘’F 

26.6 

R-500 

(  18) 

CCI2F2-CH3-CHF2 

-27“F 

23.8 

(139) 

86‘'F 

20.8 

R-502 

(  31) 

CHCIF2-CCIF2CF3 

-27“F 

24.9 

(191) 

86‘’F 

21.2 

R-717 

(  16) 

NH3 

-27“F 

39.5 

(171) 

86'»F 

37.6 

Sodium  bromide 

NaBr 

780 

43.0 

Sulfu 

ric  acid 

H2SO4 

20 

1491 

Toluene 

C7«8 

25 

68.7 

Turpentine 

Cl0«16 

25 

171 

Water 

H2O 

25 

96.1 

Zinc 

Zn 

280 

27.9 

Me 

rcury  looks  best 

at  V-1.21E-6  ft 

^/sec  at 

25®C  (77“F), 

but 

mercury  is  opaque, 

highly  reflect 

ive  and 

impossible  to 

5.  All  the  refrigerants,  R-XYZ,  are  for  saturated  liquid. 
The  first  value  is  taken  at  -27®F,  and  the  second  value  is  at 
ee^F. 
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use 


since  high  speed  photography  is  ineffective  as  are 
magnetic  switches  or  other  commonly  used  methods  of  vehicle 
monitoring.  Other  fluids  are  grouped  around  ~3.E-6  ft^/sec  at 
standard  conditions;  however,  due  primarily  to  their 
flammability  and/or  volatility  they  were  eliminated.  The 
practical  choice  was  water  at  1.07E-5  ft^/sec  at  70°F  since 
the  viscosity  is  reasonably  low  and  filming  could  be  easily 
accomplished.  Using  water  also  allows  the  use  of  aluminum 
models,  which  are  easily  produced  on  a  lathe  and  have  enough 
negative  buoyancy  to  fall  reasonably  fast  yet  are  light 
enough  to  reach  terminal  velocity  in  a  5  foot  tube. 

There  are  two  primary  concerns  when  using  water  or  any 
other  liquid.  The  first  is  cavitation  and  the  second  is  model 
stability.  Cavitation  is  undesirable  because  it  creates  a  low 
pressure  zone  (concern  is  primarily  at  the  tail  region)  that 
will  produce  a  high  profile  drag.  Cavitation  can  be  reduced 
or  eliminated  by  avoiding  flow  separation  at  the  rear  of  the 
model.  in  this  connection,  it  is  helpful  to  keep  model  speed 
down  and  simultaneously  place  a  conical  tail  on  the  model.  If 
higher  model  speeds  or  blunt  tail  cones  are  desired,  the 
problem  can  be  eliminated  by  pressurizing  the  system  so  that 
at  no  point  in  the  flow  field  does  the  fluid  experience  less 
than  the  vapor  pressure.  The  approach  taken  in  this  study  was 
to  keep  model  speed  low  (terminal  velocity  less  than  5 
ft/sec)  and  to  have  an  integral  conical  tail  with  a  slope  of 
1:3.  The  high  speed  film  was  then  closely  examined  for 
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tell-tale  vapor  bubbles  at  the  rear  of  the  falling  model, 
thus,  it  was  decided  that  cavitation  was  not  a  factor.  Model 
stability  was  provided  by  a  central  guidewire  as  discuss 

Chap.  3. 

A  presentation  of  experimental  data  of  distance  vs.  time, 
i.e.,  z  vs.  t,  from  tbe  drop  tests  is  given  in  the  following 
three  tables.  Table  A-4  presents  five  representative  tuns 
of  the  625  body  falling  through  a  water  filled  tube  of  1 
inside  diameter  (d.0.625",  D-l.O").  Table  A-5  gives  five 
runs  for  tbe  750  body  (d-0.75”,  D-l").  Table  A-6  presents 

five  runs  for  the  875  body  (d-0.875",  D-1.0"). 
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Table  A-4.  Model  625  Drop  Test  Data 


TIME 
sec . 
0.00 
0.01 
0.02 
0.03 
0.04 
0.05 
0.06 
0.07 
0.08 
0.09 
0.10 
0.11 
0.12 
0.13 
0.14 
0.15 
0.16 
0.17 
0.18 
0.19 
0.20 
0.21 
0.22 
0.23 
0.24 
0.25 
0.26 
0.27 
0.28 
0.29 
0.30 


Z(l) 

in . 

Z(2) 

in . 

Z(3) 

in . 

0.00 

0.00 

0.00 

0.03 

0.01 

0.04 

0.08 

0.04 

0.09 

0.14 

0.07 

0.17 

0.20 

0.12 

0.29 

0.26 

0.19 

0.40 

0.35 

0.32 

0.55 

0.50 

0.48 

0.67 

0.69 

0.67 

0.80 

0.89 

0.85 

0.97 

1.14 

1.01 

1.15 

1.34 

1.24 

1.34 

1.54 

1.47 

1.61 

1.74 

1.70 

1.85 

1.96 

1.93 

2.16 

2.22 

2.23 

2.42 

2.50 

2.51 

2.75 

2.80 

2.88 

3.05 

3.13 

3.18 

3.32 

3.40 

3.51 

3.68 

3.86 

3.87 

4.00 

4.28 

4.13 

4.37 

4.66 

4.54 

4.77 

5.05 

4.96 

5.19 

5.48 

5.31 

5.62 

5.91 

5.67 

6.03 

6.30 

6.09 

6.47 

6.63 

6.44 

6.83 

7.11 

6.86 

7.33 

7.52 

7.42 

7.68 

7.99 

7.87 

8.09 

Z(4) 

in . 

Z(5) 

in . 

z 

avg 
in . 

0.00 

0.00 

0.00 

0.00 

0.04 

0.02 

0.01 

0.09 

0.06 

0.06 

0.15 

0.12 

0.16 

0.21 

0.19 

0.30 

0.29 

0.29 

0.44 

0.38 

0.41 

0.59 

0.56 

0.56 

0.77 

0.75 

0.74 

0.91 

0.91 

0.91 

1.09 

1.05 

1.09 

1.32 

1.27 

1.30 

1.58 

1.58 

1.55 

1.79 

1.76 

1.77 

2.10 

2.00 

2.03 

2.33 

2.22 

2.29 

2.59 

2.59 

2.59 

2.89 

2.89 

2.90 

3.17 

3.18 

3.19 

3.48 

3.54 

3.52 

3.84 

3.95 

3.90 

4.26 

4.31 

4.27 

4.70 

4.67 

4.67 

5.13 

5.01 

5.07 

5.54 

5.38 

5.47 

5.92 

5.82 

5.87 

6.28 

6.33 

6.29 

6.72 

6.80 

6.68 

7.21 

7.22 

7.14 

7.62 

7.67 

7.58 

8.08 

8.13 

8.03 
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Table  A“4 .  continued 


TIME 

sec . 

Z(l) 

in . 

Z(2) 

in. 

Z(3) 

in . 

Z(4) 

in . 

Z(5) 

in . 

^avg 
in . 

0.31 

8.51 

8.32 

8.65 

8.59 

8.54 

8.52 

0.32 

9.01 

8.95 

9.08 

9.10 

9.09 

9.05 

0.33 

9.46 

9.45 

9.50 

9.57 

9.56 

9.51 

0.34 

9.94 

9.91 

9.93 

10.06 

10.04 

9.97 

0.35 

10.48 

10.34 

10.38 

10.57 

10.61 

10.48 

0.36 

10.99 

10.82 

10.86 

11.13 

11.11 

10.98 

0.37 

11.49 

11.31 

11.35 

11.58 

11.62 

11.47 

0.38 

12.00 

11.80 

11.89 

12.07 

12.20 

11.99 

0.39 

12.55 

12.32 

12.43 

12.53 

12.62 

12.49 

0.40 

13.10 

12.82 

12.96 

13.03 

13.13 

13.01 

0.41 

13.65 

13.32 

13.50 

13.54 

13.63 

13.53 

0.42 

14.19 

13.82 

14.04 

14.10 

14.14 

14.06 

0.43 

14.74 

14.32 

14.58 

14.67 

14.65 

14.59 

0.44 

15.30 

14.82 

15.12 

15.25 

15.17 

15.13 

0.45 

15.85 

15.33 

15.67 

15.82 

15.76 

15.69 

0.46 

16.40 

15.89 

16.23 

16.40 

16.35 

16.25 

0.47 

16.95 

16.45 

16.80 

16.98 

16.94 

16.82 

0.48 

17.51 

17.01 

17.36 

17.56 

17.53 

17.40 

0.49 

18.08 

17.56 

17.92 

18.15 

18.12 

17.97 

0.50 

18.65 

18.14 

18.49 

18.74 

18.69 

18.54 

0.51 

19.22 

18.73 

19.03 

19.32 

19.27 

19.11 

0.52 

19.78 

19.33 

19.57 

*** 

19.84 

19.63 

####  steady  state  #### 

0.53 

*** 

19.92 

20.11 

*  ★  ★ 

20.42 

20.15 

*** 

20.52 

20.65 

*** 

20.58 
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Table  A-5.  Model  750  Drop  Test  Data 

TIME  rrn  zT2i  ztt)  ztt)  ztbt 

sec.  in.  in . _ in . _ in . _ in . _ in . 


0.00 

0.01 

0.02 

0.03 

0.04 

0.05 

0.06 

0.07 

0.08 

0.09 

0.10 

0.11 

0.12 

0.13 

0.14 

0.15 

0.16 

0.17 

0.18 

0.19 

0.20 

0.21 

0.22 

0.23 

0.24 

0.25 

0.26 

0.27 

0.28 

0.29 

0.30 


0.00 

0.00 

0.00 

0.00 

0.03 

0.10 

0.20 

0.30 

0.44 

0.57 

0.71 

0.84 

0.99 

1.21 

1.48 

1.66 

1.83 
2.03 

2.27 

2.50 
2.75 
3.02 
3.26 

3.48 

3.79 
4.09 
4.35 

4.58 
4.91 
5.20 

5.46 


0.00 

0.00 

0.02 

0.06 

0.13 

0.23 

0.33 

0.40 

0.50 

0.59 

0.70 

0.86 

1.07 

1.27 
1.42 

1.60 

1.79 
2.11 
2.33 
2.53 

2.80 
3.06 

3.28 

3.61 

3.80 
4.16 
4.44 
4.70 
4.96 
5.26 
5.57 


0.00 

0.00 

0.01 

0.04 

0.08 

0.15 

0.26 

0.39 

0.53 

0.66 

0.76 

0.89 

1.04 

1.22 

1.42 
1.63 

1.84 
2.06 

2.32 

2.58 

2.84 
3.08 

3.33 

3.58 

3.84 
4.10 
4.36 

4.61 
4.87 

5.14 

5.43 


0.00 

0.03 

0.08 

0.09 

0.14 

0.23 

0.34 

0.44 

0.54 

0.65 

0.80 

0.95 

1.11 

1.31 

1.47 

1.62 
1.82 
2.07 

2.29 

2.51 

2.77 
3.03 
3.23 

3.49 
3.70 
3.96 
4.21 

4.47 
4.73 
5.01 

5.28 


0.00 

0.02 

0.06 

0.11 

0.17 

0.22 

0.33 

0.45 

0.53 

0.62 

0.76 

0.91 

1.07 

1.26 

1.46 

1.63 

1.86 

2.14 
2.32 

2.52 
2.76 
3.05 

3.28 

3.59 
3.79 
4.03 

4.29 
4.57 

4.84 

5.15 
5.45 


0.00 

0.01 

0.03 

0.06 

0.11 

0.19 

0.29 

0.39 

0.51 

0.62 

0.74 

0.89 

1.06 

1.25 

1.45 

1.63 

1.83 

2.08 

2.31 

2.53 

2.78 
3.05 

3.28 
3.55 

3.78 
4.07 
4.33 

4.59 
4.86 

5.15 
5.44 
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Table  A- 5.  continued 


TIME 

sec . 

Z(l) 

in . 

Z(2) 

in. 

Z(3) 

in . 

Z(4) 

in . 

Z(5) 

in . 

z 

avg 
in . 

0.31 

5.77 

5.83 

5.70 

5.61 

5.73 

5.73 

0.32 

6.11 

6.09 

5.96 

5.86 

6.00 

6.00 

0.33 

6.37 

6.37 

6.21 

6.12 

6.27 

6.27 

0.34 

6.67 

6.68 

6.53 

6.39 

6.55 

6.56 

0.35 

7.02 

6.96 

6.84 

6.67 

6.83 

6.86 

0.36 

7.32 

7.25 

7.13 

6.95 

7.12 

7.15 

0.37 

7.58 

7.54 

7.43 

7.23 

7.41 

7.44 

0.38 

7.89 

7.82 

7.72 

7.51 

7.70 

7.73 

0.39 

8.20 

8.12 

8.01 

7.78 

7.99 

8.02 

0.40 

8.47 

8.42 

8.30 

8.07 

8.28 

8.31 

0.41 

8.76 

8.68 

8.58 

8.38 

8.59 

8.60 

0.42 

9.06 

8.94 

8.87 

8.70 

8.90 

8.89 

0.43 

9.37 

9.19 

9.15 

9.01 

9.20 

9.19 

0.44 

9.72 

9.51 

9.44 

9.32 

9.51 

9.50 

0.45 

10.01 

9.81 

9.73 

9.61 

9.81 

9.79 

0.46 

10.27 

10.13 

10.03 

9.88 

10.11 

10.08 

0.47 

10.54 

10.46 

10.32 

10.15 

10.40 

10.37 

0.48 

10.86 

10.74 

10.62 

10.41 

10.70 

10.67 

####  steady  state  #### 

0.49 

11.15 

11.02 

10.91 

10.68 

10.99 

10.95 

0.50 

11.42 

11.29 

11.19 

10.99 

11.29 

11.23 

0.51 

11.68 

11.63 

11.46 

11.30 

11.58 

11.53 

0.52 

11.95 

11.99 

11.74 

11.62 

11.88 

11.83 

0.53 

12.23 

12.28 

12.02 

11.94 

12.17 

12.13 

0.54 

12.52 

12.57 

12.29 

12.26 

12.47 

12.42 

0.55 

12.81 

12.86 

12.57 

12.53 

12.76 

12.71 

0.56 

13.10 

13.16 

12.86 

12.80 

13.05 

12.99 

0.57 

13.38 

13.45 

13.15 

13.07 

13.33 

13.28 

0.58 

13.67 

13.74 

13.43 

13.34 

13.62 

13.56 

0.59 

13.95 

14.03 

13.72 

13.62 

13.91 

13.84 

0.60 

14.23 

14.32 

14.00 

13.90 

14.19 

14.13 
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Table  A-5.  continued 


TIME 

SBC  • 

Z(l) 

in . 

Z(2) 

in . 

Z(3) 

in . 

Z(4) 

in. 

Z(5) 

in . 

^avg 
in . 

0.61 

14.50 

14.62 

14.29 

14.22 

14.47 

14.42 

0.62 

14.78 

14.91 

14.57 

14.54 

14.74 

14.71 

0.63 

15.06 

15.20 

14.86 

14.86 

15.01 

15.00 

0.64 

15.36 

15.49 

15.15 

15.19 

15.29 

15.29 

0.65 

15.67 

15.78 

15.45 

15.48 

15.57 

15.59 

0.66 

15.98 

16.07 

15.75 

15.78 

15.87 

15.89 

0.67 

16.29 

16.35 

16.05 

16.07 

16.16 

16.19 

0.68 

16.60 

16.62 

16.35 

16.37 

16.46 

16.48 

0.69 

16.89 

16.90 

16.63 

16.66 

16.75 

16.77 

0.70 

17.19 

17.17 

16.90 

16.93 

17.04 

17.05 

0.71 

17.48 

17.44 

17.17 

17.20 

17.33 

17.33 

0.72 

17.78 

17.72 

17.44 

17.47 

17.62 

17.60 

0.73 

18.07 

18.00 

17.71 

17.74 

17.91 

17.89 

0.74 

18.36 

18.28 

17.99 

18.01 

18.20 

18.17 

0.75 

18.64 

18.57 

18.27 

18.30 

18.47 

18.45 

0.76 

18.92 

18.85 

18.56 

18.59 

18.74 

18.73 

0.77 

19.20 

19.14 

18.85 

18.89 

19.01 

19.02 

0.78 

19.48 

19.45 

19.14 

19.18 

19.28 

19.30 

0.79 

19.78 

19.75 

19.43 

19.48 

19.54 

19.60 

0.80 

20.08 

20.06 

19.72 

19.77 

19.83 

19.89 

0.81 

20.38 

20.37 

20.01 

20.06 

20.12 

20.18 

0.82 

20.68 

20.67 

20.30 

20.34 

20.41 

20.48 

0.83 

20.98 

20.98 

20.58 

20.63 

20.70 

20.78 

0.84 

21.27 

21.30 

★  *  ★ 

20.92 

20.99 

21.12 

0.85 

21.53 

21.61 

★  ★  * 

★  ★  ★ 

21.28 

21.47 

0.86 

21.79 

21.92 

*** 

★  *  ★ 

*  *  ★ 

21.85 

0.87 

22.05 

22.23 

*  *  ★ 

★  *  ★ 

22.14 

0.88 

22.31 

*** 

*  *  * 

*  ★  * 

22.31 
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Table 

A-6.  Model  875 

Drop  Test  Data 

TIME 

sec . 

2(1) 

in . 

2(2) 

in. 

2(3) 

in. 

2(4) 

in . 

2(5) 

in . 

z 

avg 
in . 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.01 

0.01 

0.02 

0.01 

0.01 

0.01 

0.01 

0.02 

0.03 

0.04 

0.03 

0.03 

0.03 

0.03 

0.03 

0.06 

0.09 

0.05 

0.08 

0.07 

0.07 

0.04 

0.11 

0.14 

0.09 

0.18 

0.15 

0.13 

0.05 

0.15 

0.19 

0.14 

0.26 

0.20 

0.19 

0.06 

0.20 

0.24 

0.24 

0.32 

0.25 

0.25 

0.07 

0.25 

0.29 

0.28 

0.42 

0.30 

0.31 

0.08 

0.34 

0.37 

0.34 

0.50 

0.35 

0.38 

0.09 

0.44 

0.47 

0.44 

0.57 

0.41 

0.46 

0.10 

0.53 

0.57 

0.53 

0.66 

0.50 

0.56 

0.11 

0.62 

0.63 

0.62 

0.73 

0.59 

0.64 

0.12 

0.70 

0.70 

0.71 

0.81 

0.66 

0.72 

0.13 

0.75 

0.80 

0.80 

0.90 

0.76 

0.80 

0.14 

0.84 

0.89 

0.84 

1.00 

0.86 

0.89 

0.15 

0.94 

0.98 

0.90 

1.09 

0.96 

0.97 

0.16 

1.03 

1.07 

0.99 

1.17 

1.05 

1.06 

0.17 

1.07 

1.14 

1.09 

1.25 

1.15 

1.14 

0.18 

1.17 

1.21 

1.17 

1.33 

1.24 

1.23 

0.19 

1.26 

1.31 

1.24 

1.42 

1.33 

1.31 

0.20 

1.35 

1.40 

1.34 

1.50 

1.41 

1.40 

0.21 

1.45 

1.49 

1.44 

1.60 

1.50 

1.50 

0.22 

1.53 

1.58 

1.54 

1.70 

1.58 

1.59 

0.23 

1.62 

1.68 

1.64 

1.79 

1.66 

1.68 

0.24 

1.71 

1.77 

1.74 

1.89 

1.75 

1.77 

0.25 

1.80 

1.86 

1.85 

1.99 

1.84 

1.87 

0.26 

1.89 

1.96 

1.95 

2.08 

1.93 

1.96 

0.27 

1.97 

2.05 

2.03 

2.18 

2.02 

2.05 

####  steady  state  #### 

0.28 

2.06 

2.14 

2.11 

2.28 

2.11 

2.14 

0.29 

2.14 

2.23 

2.20 

2.38 

2.21 

2.23 

0.30 

2.23 

2.32 

2.28 

2.47 

2.29 

2.32 
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Table  A-6 .  continued 


TIME 

sec . 

Z(l) 

in . 

Z(2) 

in. 

Z(3) 

in . 

Z(4) 

in . 

Z(5) 

in. 

z 

avg 
in . 

0.31 

2.31 

2.40 

2.37 

2.54 

2.38 

2.40 

0.32 

2.40 

2.49 

2.46 

2.62 

2.47 

2.49 

0.33 

2.48 

2.57 

2.55 

2.70 

2.56 

2. SI 

0.34 

2.55 

2.65 

2.64 

2.78 

2.65 

2.66 

0.35 

2.63 

2.73 

2.73 

2.86 

2.74 

2.74 

0.36 

2.71 

2.80 

2.83 

2.95 

2.83 

2.82 

0.37 

2.80 

2.88 

2.90 

3.04 

2.91 

2.91 

0.38 

2.89 

2.95 

2.98 

3.13 

3.00 

2.99 

0.39 

2.99 

3.03 

3.05 

3.22 

3.09 

3.08 

0.40 

3.08 

3.12 

3.13 

3.30 

3.18 

3.16 

0.41 

3.18 

3.21 

3.21 

3.38 

3.26 

3.25 

0.42 

3.27 

3.30 

3.29 

3.46 

3.35 

3.33 

0.43 

3.35 

3.40 

3.37 

3.54 

3.44 

3.42 

0.44 

3.44 

3.49 

3.46 

3.62 

3.52 

3.51 

0.45 

3.53 

3.58 

3.55 

3.70 

3.60 

3.59 

0.46 

3.62 

3.66 

3.63 

3.79 

3.68 

3.68 

0.47 

3.70 

3.75 

3.72 

3.87 

3.76 

3.76 

0.48 

3.79 

3.83 

3.81 

3.96 

3.83 

3.84 

0.49 

3.88 

3.91 

3.90 

4.04 

3.91 

3.93 

0.50 

3.97 

4.01 

3.98 

4.13 

4.01 

4.02 

0.51 

4.06 

4.12 

4.07 

4.23 

4.10 

4.12 

0.52 

4.15 

4.22 

4.15 

4.32 

4.20 

4.21 

0.53 

4.24 

4.32 

4.24 

4.42 

4.30 

4.30 

0.54 

4.33 

4.43 

4.32 

4.51 

4.39 

4.40 

0.55 

4.42 

4.49 

4.41 

4.61 

4.47 

4.48 

0.56 

4.51 

4.55 

4.49 

4.70 

4.55 

4.56 

0.57 

4.59 

4.61 

4.57 

4.79 

4.63 

4.64 

0.58 

4.68 

4.67 

4.64 

4.88 

4.71 

4.72 

0.59 

4.76 

4.72 

4.72 

4.98 

4.79 

4.79 

0.60 

4.85 

4.79 

4.79 

*  * 

*  * 

4.81 
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APPENDIX  B 


DEFORMING  AND  ACCELERATING  CONTROL  VOLUME  ANALYSIS 


A  control  volume  analysis  proceeding  from  an  inertial 
coordinate  system  requires  proper  development  of  the  integral 
form  of  the  momentum  equation  with  a  deforming  and 
accelerating  boundary,  figure  12  .  The  moving  boundary  is 
assumed  to  be  attached  to  the  nose  of  the  body  and  the 
stationary  boundary  is  fixed  to  the  plugged  end  of  the  tube. 
The  fixed  tube  wall  provides  the  remaining  boundary.  Begin 
the  analysis  with  Newton's  Second  Law  of  Motion: 


m 


a 


®  ®rel 


)  , 


where  ag  »  acceleration  of  non-inertial  coordinate  system 
(xy)  wrt  the  fixed  inertial  system  (XY), 

*  acceleration  of  a  particle  wrt  the  non-inertial 
system  (xy). 

Apply  Reynolds  Transport  Theorem  (RTT)  to,  m  ag 


m  a, 


p  dA  + 


CS 


n)  dA  , 


where  Vq  -  instantaneous  velocity  of  xy  system  wrt  XY 
system.  Note  that  Vq  is  not  a  function  of 
either  the  volume  or  the  shape  of  the  c.v., 

“  velocity  of  particles  wrt  the  c.s.  velocity. 


1.  notation:  c . v.-cv-control  volume, 

c . s . “cs-cont rol  surface,  A  •=  volume,  A  ■  area 
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-> 

m  a, 


d 

dt 


“4 

V, 


+  V, 


p  dA 
cv 

Now  use  Law  of  Conservation  of  Mass: 


p  (Vj,  •  n)  dA  . 
cs 


dt  J 


p  d  A  + 
cv 


p  (V^  •  n)  dA  , 


cs 


or 


p  (Vj.  •  n)  dA  =  - 
CS 


Expanding: 


m  a. 


dt  J 


p  dA  + 


cv 


dt  J 


dV, 


dt 


P  dA  . 
cv 


p  dA 
cv 


V, 


dt  J 


p  dA  , 
cv 


dv. 


m  ao  = 


dt 


P  dA . 

cv 

Now,  apply  Reynolds  Transport  Theorem  to  the  quantity,  m 


m  a 


rel 


dt  J 


re 


p  dA  + 


cv 


pVrei  *  ")  dA, 


cs 


where 


••  instantaneous  particle  velocity  in  xy. 


dV, 


A  Z  ?  -  1:0 

dt  J 


p  dA  + 


cv 


dt  J 


'^rel  P 
cv 


pVrei  (Vj.  •  n)  dA , 
cs 


which  is  the  law  of  conservation  of  momentum  for  a 
deforming,  accelerating  control  volume.  The  form  presented 
here  is  easily  seen  to  reduce  to  the  familiar  form  commonly 
used  in  fixed  control  volume  analysis. 
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APPENDIX  C 


EXPERIMENTAL  ACCURACY 

This  appendix  obtains  quantitative  estimates  of 
experimental  accuracy  for  the  measured  variables.  In  cases 
where  quantitative  information  is  unavailable,  qualitative 
estimates  are  made  concerning  relative  error. 

Measured  Displacement^ 

Model  displacement,  S  (or  Z),  was  calculated  from 


where  s_  »  measured  distance,  L  «  actual  model  length, 
m 

L  »  measured  model  length.  Note;  L/L„  -  conversion  scale, 
m  ^  m 

The  relative  error  is  then 

AS  As„  AL  AL„ 

m  m 


As  AL 

-  &  - -  ;  Both  terms  are  errors  introduced  by  reading  the 

^m  m 

ruler.  The  ruler  was  marked  in  .01  increments.  At  times  the 
image  of  the  model  would  become  slightly  blurred  or  out  of 
focus  when  projected  on  the  screen.  To  help  minimize  errors 
the  measurements  were  taken  twice,  by  two  different  people, 
then  averaged.  The  maximum  A  was  +.03  inches.  Since  s_  and  L_ 
were  both  in  a  range  from  3  to  4  inches,  the  pessimistic 


1.  subscript:  m  «  measured 
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estimate  of  relative  error  will  be  taken  as  +.03/3 


±1%. 


AL 

—  :  This  term  comes  from  errors  measuring  the  actual  length 
L 

of  the  model  using  a  micrometer.  The  micrometer  is  accurate 
to  at  least  +.001  inch.  A  skilled  machinist  read  the 
measurements.  The  value  for  L  ranged  from  3.5  to  7.125 
inches.  The  relative  error  is  then  +.001/3.5  =  .03%, which  is 
negligible . 

The  relative  error  is  then 

AS  As„  AL  AL„ 

m  m 

—  *  -  +  —  +  -  • 

S  L  L„ 

HI  m 

S  rms  error:  (.01^  +  .0003^  +  .01^)°’^  »  1.4%. 

Velocity  Error 

The  terminal  velocity  was  calculated  using  the  slope  of  the 
displacement  vs.  time  curve.  The  equation  for  velocity  is 


£S 

5t 

Relative  error  is 

AV  ASS  ASt 


V  5S  St 

ASS 

7^  :  This  term  may  be  described  as  "the  error  reading  SS 
65 

divided  by  SS."  The  maximum  error  reading  SS  is  2*+.02,  so 
ASS  »  +.04.  At  a  terminal  velocity  of  9  in./sec  for  the 
slowest  model,  SS=  0.1  inch  from  frame  to  frame.  But  care  was 
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taken  to  use  a  minimum  of  2  inches  for  6S  when  determining 
terminal  velocity.  The  result  is  a  relative  error  of  +2%. 


-  :  Similar  to  the  preceding  term,  this  term  can  be 

&  t 

thought  of  as  "the  error  in  reading  St."  Time  was  accurately 
measured  with  a  stroboscope  mounted  internally  to  the 
camera.  The  strobe  was  checked  against  a  laboratory  quality 
stroboscope  accurate  to  +.0001  sec,  and  found  to  have 
essentially  perfect  synchronization.  The  internal  strobe 
placed  timing  marks  on  the  film  at  6t«  .01  sec.  At  the 
highest  terminal  velocity  measured,  56  in./sec,  care  was  taken 
to  use  at  least  2  inches  of  travel,  which  takes  St-. 036  sec 
or  3.6  timing  marks.  The  most  pessimistic  approach  is  to 
assume  an  error  of  1/10  of  the  interval  between  two  timing 
marks.  The  result  is  a  relative  error  of  +.001/. 036=  2.8%.  The 
maximum  relative  error  for  terminal  velocity  is 
av  ASS  ASt 

V  5S  St 

V  rms  error;  (.02^  +  .028^)^*^  -  3.4%  . 

Drag  Coefficient 

At  terminal  velocity,  total  drag  equals  the  net  weight  of  the 
model.  Drag  coefficient  is  calculated  from  the  relation 


C 


d 


2  ^n 
p  a 
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Relative  error  is  then 


whe  re 

and 

T-  :  The  diameter  of  the  model  was  measured  using  a 
d 

micrometer.  The  maximum  relative  error  is  +.001/. 625-  +.2%. 

AD 

—  :  The  tube  was  measured  internally  to  obtain  a  diameter 
D 

at  several  places  on  both  ends.  It  was  found  that  the 
variation  in  internal  diameter  was  no  more  than  +.005  inches. 
Relative  error  is  then  +.005/1.0=  .5%. 

A3  Ad  AD 

So :  —  —  —  +  —  • 

3  d  D 

3  rms  error:  (.002^  +  .005^)^’^  -  0.54%  . 

AW 

—  :  The  models  were  weighed  on  a  laboratory  quality  balance 

w 

scale  that  was  calibrated  immediately  beforehand.  The  scale 
is  taken  to  be  accurate  to  within  at  least  +.5  gm.  The 
lightest  model  had  a  mass  of  62.47  gm.  The  relative  error  is 
then  +.5/62.5  or  +0.8%. 


^d 

AW. 

W. 


AW^  Ap 


W 


n 


AV  Aa 
+  —  +2  —  +  —  ' 
p  V  a 


n 


AW  Ap  A A 

W  p  A 


Aa 

a 


Ad 

2  — 
d 
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A  p 

—  ;  The  temperature  of  the  water  was  measured  at  the  top  of 
p 

tube  before  each  test  series.  The  accuracy  of  the  temperature 
is  +.5®F.  Density  was  then  found  from  ref  [42]  to  within  +5 
kg/m**3.  The  relative  error  is  then  +5/1000.-  .5%. 


Av  . 

—  :  Kinematic  viscosity  was  found  from  ref  [42]  to  within 

V 

a  relative  error  of  +.5%  using  the  temperature  of  the  water. 


ARe ,  AV  Ad  Av 

So:  -  “  —  +  —  +  —  • 

Re^  V  d  V 

Re^  rms  error:  (.034^  +  .002^  +.005^)®*^  -  3.4%. 


AA 

—  :  Volume  of  the  models  was  found  by  displacing  water  in  a 
A 

large  graduated  cylinder.  The  readings  were  repeated  a  number 
of  times  by  two  different  people.  The  readings  were 
consistent  to  within  ±-5  ml.  The  model  with  the  smallest 
volume  had  22.4  ml.  Relative  error  is  +.5/22.4=  +2.2%.  The 
relative  error  is  then 


AW  Ap  AV  Aa 

— -  +  —  +2  —  +  — 

p  V  a 


rms  error: 


(.024^  +  .005^  +  .068^  +  .004^)°*^ 


7.2%  , 


where 


AW  Ap  AA 

W  p  A 


W  rms  error: 
n 


(.008^  +  .005^  +  .022^)°*^  =  2.4%, 


155 


Aa 

and  —  “  2(  +  .4%)  *•  +.8% 

a 


La 


—  =  25%  estimated  from  the  stiffness  of  eq.  2-4a, 


oe 


—  =  25%  estimated  from  the  stiffness  of  eq.  2-4b. 


oe 


Relative  error; 


Am  AW 


m  W 


m  rms  error  *  W  rms  error  ■  0.8%  . 


Relative  error: 


m,. 


Relative  error: 


Ap  LA 

Ae 

—  +  —  + 

2 —  • 

P  A 

3 

(  .005^  + 

.022^  +  .01^) 

-  2.5% 

Am 

Aa 

—  +  - . — 

+  —  +  -  + 

m 

^  ^dt 

a 

a  rms  error  »(  .  008^+ . 025^+ . 005^+ . 072^+ . 004^ )  '  ^ 
o 


Ab^  AVt  Aa^ 
Relative  error;  -  »  -  +  - 


V.  a^ 
o  t  o 


b  rms  error  -  (.072^  +  .086^)  =  11.2% 
o 


Am 


Relative  error: 


he 


Ab^^  AW 
oe 

-  +  — 


■"he  "n  '^t  >=06  " 


m,  rms  error 
he 


( .024^+. 034^+. 25^+. 008^) 


=  25.4! 


8.6% 
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APPENDIX  D 


INVISCID  ANALYSIS 

This  appendix  presents  the  development  of  the  inviscid 
finite  difference  formulation  of  a  cylindrical  body  falling 
through  fluid  in  a  closed-end  tube. 

Problem  Formulation 

Assuming  an  incompressible  fluid  with  constant  transport 
properties  the  general  governing  partial  differential 
equations  are,  figure  12  : 

Continuity  (Law  of  Conservation  of  Mass,  LCM): 

^  •  V  =  0 . 

Navier-Stokes ; 

-> 

pg  -  ^p  +  7  •  T .  .  «  p  —  • 

^  Dt 

The  cylindrical  coordinates  (r,0,z)  are  related  to  the 
Cartesian  (x,y,z)  by 

X  «  r  cose,  y  =  r  sine,  z  ■  z. 

Let  V  =  Vj.1  +  Vq]  +  v^k  =  ui  +  V3  +  wk, 

where  v^  «  u  =  velocity  in  r  direction, 

=*  V  ■  velocity  in  e  direction, 

V  =  w  =  velocity  in  z  direction, 
z 

It  is  convenient  to  list  some  useful  relationships  for 
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cylindrical  coordinates : 


D 

Dt 


V  3 


7p  -  — 


V 


* 

—  +  u  - 

—  + 

w  — 

3t  3r  r 

30 

3z 

3p  ^  1 

—  1  +  — 

3r  r 

30 

iz\ 

1  3 

1  3v 

^  3w 

«: 

- (ru) 

+  —  — 

+1  — 

t 

r  3r 

r  30 

3z 

3 

1  3 

3 

« 

u  —  + 

—  V  — 

+  w 

—  / 

3r 

r  30 

3z 

3^  1 

3  1 

32 

32 

^  "7 

77^  7 

2 

30^ 

77 

Assuming  steady  axisymmetric  flow,  all  0  and  t  dependence  is 
ignored.  The  additional  inviscid  assumption  gives  *  0. 

These  assumptions  allow  one  to  write: 

Continuity  : 


1  3  3w 

—  —  (  ru )  +  —  =  0  . 

r  3r  3z 

Navier-Stokes  (ignoring  all  body  forces, i.e., gravity) : 


3u 

3u  1 

3p 

r : 

u  —  +  w 

—  «  _  _ 

3r 

3z  p 

3r 

0: 

0 

H 

0 

3w 

3w  1 

3p 

z : 

u  —  +  w 

-  n  _  _ 

3r 

3z  p 

3z 

Rewriting  continuity: 


3  3 

—  (ru)  +  —  ( rw)  =  0. 
3  r  3z 
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Introducing  the  Stokes  stream  function  \j/  with  units 


3 

3iJ/' 

3 

'dxp' 

3r 

3z. 

dZ 

.3r. 

sec 


Comparing  the  last  two  equations  it  is  seen  that 


9^  1  dll' 

ru  «  -  —  and  u  «  -  —  —  / 

3  z  r  3  z 

1  d\f/ 

rw  -  —  and  w  = - - 

3r  r  3r 


The  irrotational  assumption  is  made:  o) 


^  X  ^ 


CO 


^  X  V 


1 

1 

-f 

—  1 

D 

—  k 

r 

r. 

a 

3 

3 

3r 

39 

3z 

u 

rv 

w 

Wj.!  +  Wg 


+  w  k 

Z 


I 


which  reduces  to 


■3u 

3w' 

H 

<D 

3 

.3z 

3r. 

Substituting 


for 


u  and  w  in  terms  of  i|/(r,z) 


1 


l.It  is  of  interest-to  note  that  for  axi-symmet r ic  flow  the 
Laplace  equation  7^t|;  it  0,  although  =  0.  Further  the 

Cauchy-Riemann  relations  from  plane  flow  do  not  hold.  Plane 
flow:  V'  *  Lagrange  stream  function,  in'^/sec 
tp  =  Potential  function,  in^/sec 


3* 

3<(> 

3iJ; 

9<p 

94> 

3i|/ 

Cauchy-Riemann:  u=  —  : 

—  = 

— 

and  v«  —  : 

— 

-  — 

3x 

3x 

3y 

3y 

3y 

3x 

3 

Axi-symmetric  flow:  f  -  Stokes  stream  function,  in  /sec 

2 

^  *  Potential  function,  in  /sec 
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2 

D  \j/ 


3  \i/  1  3xJ/ 


3r' 


r  3r 


3^1)/ 


3z' 


for 


0  <  r  < 

-  -  o 

0  <  z  <  c 


with  the  body  defined  in  figure  8  :  (0  <  z.  <  L)  by; 


/(Ci^  -  (r,  -  z^,)2) 


r . 
1 


from  0  < 

from  r.  <  z.  <  (L  -  L,) 


r^  “  r^(L  —  from  (L  ~  ^3^  —  ^b  —  ^  ^ 

subject  to  the  following  (mixed)  boundary  conditions: 
a.)  Dirichlet:  -  0  axis  of  symmetry  (AOS)  and  on  the  body.' 


b.)  Neumann:  v. 


w 


specified  at  inlet. 


c. )  Dirichlet:  V/  =  constant  (^0)  on  the  tube  wall. 

d.  )  Neumann:  v  »  w  »•  specified  at  outlet. 

Each  of  these  boundary  conditions  are  now  discussed: 

Let  «  body  velocity,  in. /sec. 

a. )  -  0  and  '('(  Zj^)  *  0. 

b.  )  Assuming  the  inlet  velocity  is  parallel  to  the  AOS: 

2 


'Kr,0) 


^  r2 


for 


(0  <  r  <  r^)  . 


(Footnote  1  continued  from  previous  page) 


3  + 

34> 

1 

3\j; 

34) 

34> 

1 

9\|/ 

u«  —  : 

—  « 

-  — 

—  and 

V  “  w=  —  ; 

—  = 

— 

3r 

3r 

r 

dz 

^  3z 

3z 

r 

3r 

Lines  of  constant  V/  are  orthogonal  to  lines  of  constant  ^ 
in  plane  and  axi-symmetric  flow. 

2.  Mathematics  requires  that  \fi  -  constant  on  a  streamline. 
It  is  simply  convenient  to  choose  the  constant  to  be  zero 
since  the  axis  of  symmetry  and  the  body  itself  form  what 
turns  out  as  a  streamline. 

3.  tj;  must  equal  a  constant  along  the  tube  wall  because  the 
tube  is  a  streamline.  However  jp  must  be  different  from  zero 
or  else  there  would  be  no  flow  in  the  region  of  interest. 

4.  To  specify  that  v^,  =  u  »  0  malces  ^(r,z)  =  f(r)  only. 


160 


c. )  r  2  for  (0  <  z  <  c). 

2 

d. )  Assuming  the  outlet  velocity  is  parallel  to  the  AOS: 

t*/(r,c)  «  r  for  (0  <  r  <  r^). 

2 

In  the  event  that  the  outlet  is  chosen  to  be  at  some 
point  on  the  body  where  u  *  0/  say  in  the  middle  of  the 
horizontal  section,  at  z  *  d: 

V  r2  -  r.2 

4-(r,d)  -  —  r  2  <  r  <  r  ).  ; 

2  -  r/ 

O  1 

The  result  is  a  transformation  from  the  original  mixed 
boundary  value  problem  to  a  Dirichlet,  or  first,  boundary 
value  problem. 

The  next  section  will  discuss  the  discretization  and 
methods  used  to  solve  the  inviscid  relations  just  derived. 
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V 


Discretization  Method 

For  convenience  the  differential  equation  is  repeated  here; 


2  ^ 
D^«J/  =  - =• - 

—  +  — 2  “  ^ 

for 

0 

<■  r 

<  r 

ar^  r 

3r  3z 

0 

<  c 

The  terms  of  which  may  be  expressed  in  a  direct  analogous 
finite  difference  form  ^ (refer  to  figure  19): 


a^t|/ 

- 

'*1,1  * 

i|/  .  .  . 

^1 / 3-1 

ar2  ■ 

f 

3^4' 

- 

2*1, j  * 

3z2  ■ 

Az2 

1  3i|/ 

*1,1+1  ■ 

"  ^,1-1 

r  3r 

2r. 

1 

.Ar 

2 

After  substituting  into  D  f  ^  0  and  rearranging,  a  working 
formula  is  found; 


♦i,j  -  *  '♦ifi.j  *1-1, j'  *  ^  ' 

*  ”  ar)a2^ 

where  A  =  - = - =—  ,  B  =  - ^  ^ 


and 


2  2  ' 
2(Ar^  +  Az"^) 


(2r.  .  +  Ar)Az 

C  =  - - 

4r.  .(Ar^  +  Az^) 

^  f  J 


4r.  .(Ar^  +  Az^) 

^  /  J 


Residual;  R.  ^  R - 4  -  4/^  i  ^ )  -  4ii/.  .  , 

lO  2j 


where 


^1,3+1  1,3-1  1+1, J  ^1-1,3  • 


5.  i  =  column  =1,2,3,  ....  ,  MMAX 
j  =  row  *0,1, 2, 3,  ...  ,  NMAX 

Ar  =  uniform  grid  size  in  r  direction,  inches 
Az  =  uniform  grid  size  in  z  direction,  inches 
4  =  j*Ar  =  radial  distance  of  mesh  point  from  AOS 

1  /  J 
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Implementation  of  the  boundary  conditions  in  the  discretized 
form  is  obvious  and  will  not  be  discussed.  The  reader  is 
referred  to  Appendix  E.2,  showing  the  FORTRAN  implementation 
of  the  Gauss-Seidel  iterative  solution,  in  program  PSI.FOR. 
The  overall  flow  chart  used  in  PSI.FOR  is  shown  below. 


Initialize 


Generate  Boundary 
Values 


+ 


Calculate  Values 
of  ^  on  interior 
points 


- - V- - 

Calculate  Residuals 
on  interior  points 


- V - 

Maximum  Residual 


> 


- V - 

Write  out  converged  field 


Max  Res 
Greatei 
Preset 

- i 

sidual 
:  than 
Value 
\ - 1 

An  alternative  solution  may  be  achieved  by  use  of  the  TDMA  or 
Tri-Diagonal-Matrix  Algorithm  from  [37].  The  first  step  is  to 
write  the  equation  to  be  solved  in  the  proper  format;^ 


z  direction: 

or  simply  : 
( j=constant ) 


a .  «i/. 
ai  V;. 


‘^i  j  ^i  '^i-l,j  "^i  ' 

bi  +  c.  +  d.  , 


Vb  2 

where  it  is  given  that  4/.  =  -  ,  on  the  inlet,  is  known 

^  2 


V 


and 


MMAX 


l_  ^ 

—  r^  ,  on  the  outlet,  is  known. 


6. 


Subscripting  done  as  in  previous  discussion. 


Then, 


\l>.  -P.  .  +  Q., 

1  ^1+1  1  ' 


where 


b. 

1 


-  'i  ^-1 


and 


<^1  '^1  °i-i 

®i  ■  'i  ^i-1 


The  method  is  to  forward  calculate  Pj^/Q^  then  back 
calculate  along  a  horizontal  line.  Because  d^  contains 
contributions  from  above  and  below  this  line,  it  is  wise  to 
repeat  the  process  in  the  r  direction  where  contributions 
from  the  right  and  left  of  each  node  come  in  through  d^. 
Convergence  is  speeded  up  by  this  alternate  sweeping 
technique.  An  additional  aid  to  fast  convergence  is  the 
introduction  of  a  relaxation  parameter,  a.^ 


*  * 

’I'j  j  j  -  —  Rj  j  /  where  the  denotes  the  update. 

^»J  ^  r  J  4 

The  author  found  that  TDMA  with  over-relaxation  (a  =  1.2) 
converged  3  times  faster  than  the  point  by  point  Gauss-Seidel 
method.  So  as  not  to  burden  the  reader  with  the  tedious 
algebraic  details  associated  with  TDMA  a  copy  of  the  TDMA 
subroutine,  called  TDMA. SUB,  written  in  FORTRAN,  is  included 
in  Appendix  E.3.  Should  the  reader  wish  to  use  TDMA. SUB,  it 
is  simply  necessary  to  replace  CALL  PSI  with  CALL  TDMA  on 
line  #  132  in  program  PSI. FOR.  TDMA  is  carried  out  in  an 
alternate  sweeping  pattern.  An  extra  pass  through  the  annular 
region  aids  convergence  by  carrying  information  to  the  center 
of  the  flow  region  as  quickly  as  possible.  The  overall  flow 


7.  A  change  of  4  units  in  the  residual  produces  a  change  of 

-1  units  in  the  value  of  4;:i.e.,SR.  .  =  -48\j;. 

1/3  1/3 
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chart  for  TDMA  is  the  same  as 
previously  given. 


for  Gauss-Seidel,  which  was 


Hydrodynamic  Mass 

Once  the  field  is  converged  it  is  routine  matter  to 

determine  the  velocity  and  pressure  fields.  The  actual 

computations  were  carried  out  in  program  RESULTS. FOR  in 

Appendix  E.4.  Hydrodynamic  mass,  m^,  is  defined  for  constant 

0 

acceleration  in  an  inviscid  fluid,  by 


1 

2 


m. 


(w^-  V^)  +  u^ 


SA  , 


all 

fluid 

elements 

where  5A  *=  2it  r  Ar  Az,  A  *=  volume. 

An  elementary  assumption  of  constant  acceleration 
inviscid  fluid  happens  to  give  a  good  estimate  of  mj^ . 


in 


the 


Continuity ; 

Then, 

Example; 


m  =  p  (A-a)  L  ,  where  L  =  —  • 
g  e  e  g 

(A-a)  Vg  -  a  Vjj. 


'  a  " 

^  1 

<A“a> 

n2  a'2- 

-  d  J 

9  (B^  - 

1.6,  Lg»  5.33  in..  A*  .7854  in^,  a=  .3068  in.^ 

slugs 


Then,  m  =  0.00286  slugs,  p  =  0.001123  - 5 — ' 

5  in^ 

■  0.00118  slugs  (elementary  formula). 


8.  Subscripts:  g  »  gap  (between  body  and  tube),  b  ■  body 
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irij^  *  0.00121  slugs  {I  formula  from  program  RESULTS). 
Thus, the  elementary  formula  is  only  about  2.5%  low. 


Determination  of 
- P 

The  preceding  sections  have  described  the  process 
utilized  to  determine  the  tf/  field.  Then,  the  velocity  field 
was  found  using  the  definitions  of  u  and  w  in  terms  of  ij/. 
Hydrodynamic  mass  was  calculated  next,  using  u  and  w.  The 
only  remaining  flow  field  to  be  determined  is  pressure,  which 
can  be  done  through  the  steady  state  Bernoulli  equation. 
Start  with  a  streamline  drawn  from  the  inlet  to  some  point 
downstream,  figure  20;  neglect  potential  energy  changes: 


Ap 

0.5pVj^^ 


by  definition 

(Vj^  =  abs.  vel.  of  body) 


and 


—  +  —  where 
P  2 


u 


+  w 


Then,  c 


1  - 


'b  " 


The  flow  field  surrounding  body  625  was  analyzed  in  program 

PSI.FOR,  Appendix  E.2,  and  program  RESULTS. FOR,  Appendix  E.4, 

was  used  to  produce  the  velocity  field.  It  is  then  possible 

to  calculate  C  from  the  above  relation.  Figure  21  is  a 
P 

plot  Cp  for  body  625,  which  is  independent  of  Vj^  and  hence 
Reynolds  number. 
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APPENDIX  E 


COMPUTER  PROGRAMS 


E,1  UNSTEADY  FORTRAN  PROGRAM 


IMPLICIT  REAL  (A-H,J-Z) 

INTEGER  M,NUM,K 
DIMENSION  F(2) , Y(2) 

CHARACTER*80  CHOICE 

OPEN (1, NAME- 'T. DAT' ,TYPE-'NEW' ) 

OPEN (2, NAME- 'D. DAT' .TYPE- 'NEW' ) 

OPEN( 3, NAME- 'CD. DAT' ,TYPE-'NEW' ) 

M  -  0 

TIME  -  0.0 
TIMEMAX  -  1.0 
STEP  -  0.01 
NUM  -  2 
IPRINT  -  1 
I  STOP  -  0 

G  -  32.2  I  FT/SEC/SEC 

RHO  -  1.94  I  SLUGS/FT* *3 

NU  -  1.07E-05  1FT»«2/SEC 
C  FOR  625  MODEL: 

DRYWT  -  0.1640  ILBF 

FNET  -  0.105  ILBF 

DM  -  0.625  lINCHES 

L  -  5.0  lINCHES 

D  -  1.0  lINCHES 

HYMASS  -  0.00118  ISLUGS 

MASS  -  DRYWT/G  ISLUGS 

TOTMASS  -  MASS  +  HYMASS 
DM  -  DM/12. 

L  -  L/12. 

D  -  D/12. 

B  -  D/DM 
DE  -  D  -  DM 
A  -  .7854*DM»DM 
GAM  -  (B*B  -  1. )/B/B 
FUNl  -  ( (B*B+l)/2/(B*B-l) )**2 
C  Y(l)  IS  THE  INITIAL  VELOCITY  FT/SEC,  AT  TIME-0. 

C  Y(2)  IS  THE  INITIAL  POSTIION  FT,  AT  TIME-0. 

Y(l)  -  0.001  1 0.0  CAUSES  PROBLEMS 

Y{2)  -  0.0 

C  BEGIN  THE  MAIN  LOOP 

8  IF(P0S*12.  .GE.  20.0)THEN 

ISTOP  -  1 
GOTO  21 
END  IF 

6  CALL  RUNGE{NUM,Y,r, TIME, STEP, M,K) 

GOTO(10,20),K 
10  CONTINUE 

V  -  Y(l) 

POS  -  Y(2) 

0  -  .5*RH0*V»V 

RED  -  V*DM/NU 

REW  -  B*B»RED/{B+1) 

FUN2  -  .001  +  2.8*(ALOG10(REW) )**-3.1 

FUNS  -  FUN1*FUN2 

CHI  -  1.477  +  .2195*ALOG10{RED) 

ZETA  -  1.005*B**CHI 

TOS  -  FUNS*(B+1 )/4/(B+ZETA) 

TIS  -  ZETA^TOS 
TI  -  TIS*Q 
DP  -  FUNS*Q/DE 

FOR  -  3.14159*DM»L*TI  +  A»L»DP 
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s 


CS  -  FOR/A/Q 

CD  -  (GAM*GAM+(1-GAI1)**2)/GAM/GAM  +  CS 
DRAG  -  A*CD*Q 

ACC  -  (FNET  -  DRAG)/TOTMASS 
F(l)  -  ACC 
F(2)  -  V 

IF(TIME.EQ.O.O)THEN 
C  •  SET  UP  OUTPUT  TABLE  HEADINGS 

TYPE  TIMECSEC)','  VELOCITY  POSTION  CD' 

WRITE(5,23)TIHE,V»12. ,POS*12. , CD, ACC 
WRITE(1,27)TIME 
WRITE(2,27)POS»12. 

WRITE( 3 , 27 )CD 
END  IF 
GOTO  6 

C  END  MAIN  LOOP 
C  THIS  IS  THE  OUTPUT  SECTION 

20  ICOUNT  -  ICOUNT+1 

IF( INT( ICOUNT/IPRINT) ‘IPRINT  .EQ.  ICOUNT)THEN 

21  WRITE( 5,23)TIME,V»12. ,POS»12. ,CD,ACC 
WRITE(1,27)TIME 
WRITE(2,27)P0S*12. 

WRITE(3,27)CD 

23  F0RHAT(1X,5(F9.3,1X) ,2(F13.1,1X) ) 

27  F0RMAT(1X,E12.6, 'S' ) 

END  IF 

IFdSTOP  .EQ.  DSTOP 
C  GO  BACK  TO  THE  TOP  OF  THE  MAIN  LOOP 
GOTO  8 
END 

C  *****************»*»***t»»*ti,t*»t»tttt,i,»tt»ti,tt,tti,t»*» 

SUBROUTINE  RUNGE ( N , Y, F , X , H , M , K ) 

C  THIS  ROUTINE  PERFORMS  RUNGE-KUTTA  CALCULATIONS 
C  BY  GILLS  METHOD  (FROM  PROF.  FRANK  WHITE'S  ’’VISCOUS 
C  FLUID  FLOW*  BOOK). 

DIMENSION  Y(10),  F(10),  Q(10) 

M  -  M+1 

GO  TO  (1,4,5,3,7),M 

1  DO  2  I  -  1,N 

2  Q(I)  -  0. 

A  -  .5 
GO  TO  9 

3  A  -  1.70710678 

4  X  -  X+.5»H 

5  DO  6  I  -  1,N 

Y(I)  -  Y(I )+A*(F(I)*H-Q(  I) ) 

6  Q(I)  -  2.*A*H*F(1)+(1.-3*A)*Q(I) 

A  -  .2928932188 

GO  TO  9 

7  DO  8  I  -  1,N 

8  Y(I)  -  Y(I)+H*F(I)/6.-Q(I)/3. 

M  -  0 

K  -  2 

GO  TO  10 

9  K  -  1 

10  RETURN 

END 
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36 

37 

38 
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E.2  INVISCID  FORTRAN  PROGRAM 

Fortran  Program  PSI.FOR 

C  SOLVE  AXISYMMCT8IC  POTENTIAL  PLOW  AROUND  MODEL  IN  A  TUBE 
C  I-COLUHN  J-ROH  P(I,J)«VALUE  Of  PSI  e.g.  STORES  STREAM  FUNCTION, 

C  Nritttn  by:  Richard  f.  Hubbell  (July  1987)  at  Univ,  of  R.I. 

C  For  partial  coaplction  of  PhD  Mechanical  Engineering 
C  Major  Prof^  Dr.  frank  it.  white 
LOGICAL  FLOW 

CHARACTER*20  VIEW,VIEKRES 
DIMENSION  IBOOy(101),IBODY2(10n 
COMMON  /GRID/JBOO,NMAX,MMAJC 

COMMON  /RES  I  DUAL/RES  ,  IMAT  ,  JMAX  ,  RESMAX  ,  DELMAX  ,  I DELMAX ,  JDELMAX .  RELAX 
COMMON  /CENERAL/DR,OR2,D2,DZ2.I.  J.VIN,VGAP,RA,RB 

COMMON  /MATRIX/AUOl)  .8(101)  ,C(  1 0 1  )  ,  R ( 1 01 )  ,  P (  1 2 50 , 1 0 1 )  ,  FLOWC  125L*^ 

DATA  ISTAGl ,  ISTAG2  , NMAX , MMAX/401 ,684,101. 1250/ 

DATA  08,02/0.005,0.025/ 

DATA  VlN/400.0/ 

DATA  P,A,B,C/126250*0.,101*0.,101*0. ,X01*0./ 
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C 

c 


9 


OPEN(  1 , TYPE- ' NEW" . NAME- ' RES , OUT ' ) 

TYPE  ‘.'READ  IN  PSI  FROM  A  FILE?  Y  OR  N' 

ACCEPT  60, VIEW 
.F(VIEW  .EQ.  'Y')THEN 
OPEN( 2,TYPe-'OLO' . NAME- ' TOMA. PS I ' ) 

READ(  2.  •  )  I  COUNT.  ICOUNTRESMAX,  IMAX,  JMAX .  RESMAX 
READ(2,83)(  (P(I.J)  ,I-l,MMAX)  ,J-l,NMAX) 
FORMAT(12625(10(F11.8,1X) ./) ) 

TYPE  INPUT  FROM  TDMA.PSI' 

TYPE  65,((P(ISTAGWR,4  +  L),  K— 3,3),  L-3,-3.-l) 
TYPE  65,((P(ISTAG2-*.R,4>L),  R— 3.3).  L-3.-3,-l) 
TYPE  389,399' 

TYPE  9,(  iP(K,L),K-389.399) ,L-4,2,-l) 

TYPE  400,410' 

TYPE  9, ( {  P(  K,L)  ,8-400,4 10) ,L-4. 2,-1 ) 

TYPE  680,690' 

TYPE  9,(  (P(R,L),K-6a0,6  90),L-4,2.-l) 

TYPE  •, '  691,701' 

TYPE  9,( (P(K,L),R-691.701),L-4.:,-l) 
rORMAT( 3( lX.liF7.4,/) ) 

ELSE 

VIEW  -  'N' 

END  IF 


RA  •  .625/2  IRA  -  BODY  RADIUS  INCHES 

DR2  •  DR*DR 

D22  •  DZ-DZ 

RB  -  (NKAX-l)*DR 

JBOD  -  INT(RA/DR)  ♦  1 

RA  -  JBOD-DR 

RB2  •  RB*RB 

RA2  •  RA*RA 

ZSTAGl  -  ISTAGl*OZ 

2TAIL  -  1.500  'TAIL  tS  1.5  INCHES  LONG 

2STAG2  -  ISTAG2*OZ 

VGA?  -  VIH*RB2/(RB2>RA2) 


DO  52  J  -  l.NMAX 
DO  52  I  -  l.MMAX 
52  FLOWd.J)  -  .TRUE. 

C  FIRST  DEFINE  THE  NOSE  FOR  IBOOY(  }  THEN  DEFINE  THE  TAIL  FOR  IBODYZC) 


DO  50  J  -  JBOD,2,-l 

ZNOSE  •  SQRT(.RA*RA  -  (  (  J-1  J  *  DR  )  •  ♦  2  ) 
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64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 
.  76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 
109 

no 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 


IBODY(J)  •  rNT{(2STAGl  ♦  RA  •  2NOSE)/DZ) 

ZSLOPC  -  ( J-1 ) •OR*ZTAIL/RA 
50  IBODY2(J)  -  INT((2STAG2  -  2SLOPE)/D2) 

DO  53  J  ■  2,JBOO 

DO  53  I  •  IBOOY( J) ,lBOOY2( J) 

53  fLOH(nJ)  -  .fALSE. 

C  PCRrORrt  ITERATIVE  SOLUTION 
R{1)  -  0. 

DO  4  J  •  2,N?1AX 

R(J|  .  iJ-n*0R 

TOPS  «  (2.*R(J)  -  DR)*D22 

TOPC  -  (2.»RtJ)  4.  DR)*D22 

DEN  -  4.*R(J)«(D22  ♦  DR2 ) 

A(J)  -  DR2/2./( DZ2>DR2) 

B(J)  •  TOPB/DEN 
C(J)  -  TOPC/DEN 
4  CONTINUE 

tr( VIEW  .CO.  'N' )THEN 
C  INITIAL  AND  BOUNDARY  CONDITIONS 
DO  1  I  •  1  ,  HMAX 
P(I,NnAJC)  -  0.5*V1N*RB2 
1  p(r,X)  -  0.0 

^  2  J  -  JBOO^l  ,NMAX-1  !R(  J)-(  J-1  )  •OR  ;  RA-JBOD-OR  ;  SO  ADD  1  T®  ^SOD 
DO  2  I  •  IBOOY( JBOD) , IBOOY2( JBOO) 

^  -  RA2)/(RB2-RA2)  !VIH  HELPS  IT  CCf^VCAf-f 

DO  3  J  •  l.HMAX 

P(1,J)  -  0.3*VIN*R(J)*R(J) 

3  Ptnnxx.j)  .  p(i.j) 

ICOUNT  -  0 
ICOUNTRESKAJI  -  0 
ENDir 

TYPE  INPUT  IRES  (CONTROLS  EREQUENCY  OP  OUTPUT)’ 

ACCEPT  ‘.IRES 
70  CONTINUE 

TYPE  •,'VHAT  RELAXATION  PACTOR  (0  TO  3)7' 

ACCEPT  •.RELAX 

IP'RELAX  .LT.  0.  .OR.  RELAX  .OT.  3 .  I  GOTO  70 
TYPE  •,'HCW  .'RANY  (BORE)  ITERATIONS?' 

ACCEPT  ‘.ITERfLAX 

IP(  ITERftAX  .LT.  1  )THEN 

TYPE.*,'  OUTPUT  TO  PSI.PSI' 

C  TYPE  65.  (  (P(  ISTAGUR.<*L).  R— 3,3),  L-3,-3.-l) 

C  TYPE  65,  (  (P(  ISTAG2*IC,4*L),  K— 3,3),  L-3,-3.-I) 

TYPE  *,'  389,399’ 

TYPE  1.( (?(R,L  ),. 1-389, 399), L-4, 2,-1) 

TYPE  *.'  400,410' 

TYPE  9. { ( P( R,L) ,R-400, 410),L-4.2.-l) 

TYPE  *,'  680,690' 

TYPE  9,( (P(K,L)  ,  K-680, 690), L-4, 2,-1) 

TYPE  *.'  691,701’ 

TYPE  9,( (P(R,L  ),R-691, 701), L-4, 2,-1) 

TYPE  ♦,'OUTP'JT  TO  A  PILE?  Y  OR  N' 

ACCEPT  60, VIEW 
N  IP(VIEW  .EO.  ’Y')THEN 

OPEN( 3 .TYPE- 'NEW' ,NAHE-'PSI .  PSI '  ) 

WRITE! 3, • ) ICOUNT, I COUNTRESBAX, IBAX , JKAX , RESBAX 
WRITE(  3, 81  )  (  ( P(  I , J)  ,  I-1,BBAX) , J-l,NBAX) 

81  rORBAK IX, 12625 ( 10 ( Pll . 8 , IX ) ,/ ) ) 

E.NDIP 
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134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 
146 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 
167 
166 

169 

170 

171 

172 

173 

174 

175 
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178 
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182 
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184 

185 
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186 
189 


STOP  'END  or  RUN' 
ENDIf 


00  7  iter  •  l.ITERnAX 
DELrtAX  -0.0 
CALL  PSX 

C  TYPE  65,  (  |P(  ISTAGWK.  4-i.L)  , 
C  TYPE  55,( (P(ISTAC2>K,4^L), 
65  rORHAT(7(lX,7Pll.9./) ) 


K— 5,1),  L-3,.3,-1) 
K— 5,1),  L-3,.3.-l) 


62 

63 


ICOUNT  -  ICOUNT  ♦  1 

TYPE  62,  ICOUNT,  ITER,  I TERNAX  ,  ICOUNTRESNAJC ,  RESftAJC ,  IftAX,  JNAX 
TYPE  6  3,DELrtAX,lDELAAX,JOELnAJC,  RELAX 
ir(  INK  ICOUNT/IRES ) •IRES  .EQ.  ICOUNT)THEN 
RESAAX  •  0 


CALL  RESIDUALS 
ICOUNTRE3NAX  -  ICOUNT 

TYPE  62, ICOUNT, ITER, ITtRflAX,  ICOUNTRESKAX , RESRAX ,  IKAX, JRAX 

ICOUNTRES«AX,RESnAX,inAX,JnAX 
WRITE  (  l,63)DELnAX,  IDELrtAX,  JOELKAX 
rORnAT(lX,4( 14, 3X|,ri2,4,3X,2( I4.3X) ) 
rORRAT(lX,E12.5,3X,2(I4,3X),r5.J) 

END  IP 


ir(  ICOUNT  .EQ,  DTHEN 

TYPE  •,'SPOT  CHECX  BOUNDARY  CONDITIONS?  Y  OR  N' 
ACCEPT  60, VIEW 
ENDIf 

VIEW  -  'N' 

If( VIEW  .EQ.  'Y' )THEN 
TYPE  •,* INPUT  IVIEW.JVIEH* 

ACCEPT  •,IV1EW, JVIEW 
ENDIf 


iriVIEW  .EO.  'Y')THEN 
DO  31  I  .  l,nf«AX 
in  :NT(  I/IVIEH)*IV1EW.EQ.  DTHEN 
TYPE  •,(I,(PII,L),  L-l,3)) 

TYPE  * , I  I , ( P( I,L) ,  L-NnAX-2,NMAX)  ) 
TYPE  •,'  ' 

ENDIf 

31  CONTINUE 
TYPE  103 

103  rORftAK///) 


DO  32  J  -  DNRAX 
If ( INT(J/ JVIEW) •JVIEW. EQ. J)THEN 
TYPE  •,^J,(P(L,J),  L-l,3)  ) 

TYPE  •,{J,(P(L,J),  L-NnAX-2 ,RKAX)  ) 

ENDIf 

32  CONTINUE 
ENDIf 

7  CONTI.NUE 
60  rORAAKAlO) 

TYPE  62, ICOUNT, ITER,  ITERrtAX , I COUNTRESNAX , RESNAX , I?IAX , JNAX 
WRITC(  1,62)  ICOUNT,  ITER,  ITERNAX,  ICOUNTRESRAX ,  RESRAX ,  IflAX ,  JRAX 
WRITE(  I  ,  6  3  )DELRAX,  IDELRAX,  JDCL?IAX,  RELAX 
CALL  VELOCITY 

GOTO  70 
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V 


190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
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204 
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206 

207 

208 

209 

210 
211 
212 

213 

214 

215 

216 

217 

218 

219 

220 
221 
222 

223 

224 

225 

226 
::7 
:28 

229 

230 

231 

232 

233 

234 

235 

236 

237 
2  38 

239 

240 

241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 
261 
262 

263 

264 

265 


2.«P(I,J)  ♦  )/DR2 

2.*P(I,J)  ♦  P( J/D22 
I, J-1 ) )/2./R(J)/DR 


ZHD 

SUBROUTINE  RESIDUALS 
logical  PLOH 

COHMOM  /CB  1 .  D  J  .  DI 2  ,  I .  J  ,  VIH,  VGAP  ,  RA,  RB 
COnnON  /GRIO/JBOD,  NPLAX,  flKAJC 

CoJwON  yJlATRIlyAf^fnf  >  •  ^^SKAA ,  OELHAA ,  I DCLflAX  ,  JOCLRAX,  BELAID 

COMAON  /"ATRIVA(l01),B(101),C(10n,*(101),P(1250.l01).rLOH(12SO,  JCi; 
DO  I  r  • 

DO  1  J  -  2,NnAX-l 

ir(rLOK(i,j) )then 

OELSOR  •  (P(1,J*1) 

DELSOZ  -  (P(I*1,J) 

READER  -  (P1I.J*U  _  . . . 

RES  •  DELSQR  ♦  DELSQZ  -  RENDER 

t  R^SNAJC)  )THEN 

XnAA  •  I 
JNAJC  •  J 
RESKAX  -  RES 
ENDCr 

ENDir  ' 

1  CONTINUE 

TYPE  6  5  ,  (  (  P  (  INAJC>  K  ,  JKAJC^  L  )  ,  K«-3  3)  L«3  -3  -1) 

»^RITE(l.65)((P(inAX.R.jnAX;LK  r::  ' 3)  L^]'  ,, 

end 

SUBROUTINE  VELOCITY 
LOGICAL  FLOW 

CONKON  /GCNERAL/DR,  DR2  ,  D2 . OZ2 , I , J. VIN  VGAP  RA  RR 
COfWON  /GRID/JBOD,NnAX^NNJlx^^^^^ 

DENSITY  LBA-SEC •  •  2/ IN*  •  < 

DO  I  I  -  2,AAAX-1 
ir(rt,OH(I,J)|THEN 

SM*rr  rr'***  ^  ^  ^  > /2  . /R  (  J  I /DZ 

Bic«  f  ;i^^;’;'j;i|5.R(J).OR.DZ*ROW.,VZAflS*VZABS  ♦  VR-VR) 

CNDir 

1  CONTINUE 

BODYNASS  -  0.105/336.1 
HYNASS  •  BIGKE/VIN/VIN 
TYPE  10  ,flIGKE  ,  3Y?4ASS  ,  BOOYPIASS 
>fRITE(  1 ,  10  )  0IGRE,  HYXASS  ,  flODYNASS 

END 

SUBROUTINE  PSI 
LOGICAL  PLOK 

COAAON  /general/dr. DR2,DZ,DZ2, I, J.VIN, VGAP  RA  RB 
CCAMGN  /CRID/JBOO.NAAJt.AAJUC  .  VGAP ,  RA .  RB 

COAAON  /"ATRIx/A(10I),B(I01.'.C(101),R{101),p(  1250  101)  r'rwdjc/)  let) 

«AAON  ^^3I^t.RES,IAAX..«AX.RESAAX.DE;A;x':ioE?Ai5!joEU^ 

J  4pNNAX-l  lOVERALL  ITCRATICN 


IDON'T  ITERATE  THE  BOUNDARIES 


00  1  t  -  2..NMAX-1 
ir(rLOW(l,J) )TH£N 

ir(ABS(OEL)  .CT,  AflS(DELKAJC)  )THEN  ' ^  J 

DELNAJ  -  DEL 
IDELftAX  •  I 
JDELNAX  •  J 

ENDir 

P(IpJ)  -  P(I,J)  ♦  RELAX*OEL 
ENDir 
CONTINUE 
RETURN 
END 
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E.3  TRI-DIAGONAL  MATRIX  ALGORITHM  FORTRAN  PROGRAM 
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C  Written  by  Richard  r.  UubbeU  lJuly  1987  )  Univ,  ot  R.I. 

SUBROUTINE  TDflA 
LOGICAL  fLOW 

C  PCRfORrt  TOnA  IN  Z  DIRECTION  (BUT  ONLT  IN  THE  ANNULUS) 

COMMON  /GRID/JBOO,NMAJt,MMAA 

COMMON  /GENERAL/DR ,  DR2  .  D2  ,  DZ2  ,  I ,  J  ,  VIN ,  VGA? ,  RA ,  RB 

COMMON  /MATRIX/A(  101  )»B(101).C(101),R(101),P(1250,101),  PLOWl  125<?  /#/) 
COMMON  /RESIOUAL/RES,  IMAX  ,  JMAX ,  RESKAJC  ,  DELMAJC ,  X  OELMAJC  ,  JDELMAA ,  RE  li  X 
DIMENSION  P2(1250), 01(1250)  ^ 

C  BEGIN  VERTICAL  TDMA  (DO  IT  ONCE) 

DO  40  L  •  1,1 
DO  40  I  -  2,MMAX-1 
PZ(1)  -  0.0 
OZ(l)  -  P(I,1) 

DO  37  J  •  2,NMAX-1 
ir( ,NOT.PLOW{ I,J) )THEN 
ADUM  -  l.E-15 
BDUM  •  l.E-15 
CDUM  •  l.E-15 
ELSE 

ADUM  -  A(J) 

BDUM  -  B(J) 

CDUM  -  C(J) 

ENDir 

DCNOM  -  1.  -  CDUM*P2(J-l) 

P2(J)  •  BDUM/(1.  -  CDUM*PZ( J-1  )  ) 

OZNUM  -  ADUMMPlI^l.J)  ♦  P(I-1,J)) 

37  0Z(J)  -  (OZNUM  ♦  CDUM*QZ( J-l) )/DENOn 
DO  38  2  -  NMAX-1,2,-1 

DEL  -  ?2{ J) •?( I ,J+1 )  f  0Z{J)  -  P(I,J) 
ir(ABS^DEL)  .GT.  ABS  {  DELMAJt )  )  THEN 
DELMAX  •  DEL 
IDELMAX  -  I 
JDELMAX  •  J 
ENDir 

38  P(I,J)  •  P(I,J)  4.  RELAX»DEL 
40  CONTINUE 

^  Z  (BEGIN  HORIZONTAL  TDMA  (DO  IT 

PI (  I )  -  0.0 

OZ(  1)  -  P( 1. J) 

DO  7  I  •  2.MMAX-1 

ir( .NOT.rLow( i,j) )then 
ADUM  -  l.E-15 
BDUM  -  l.E-15 
CDUM  •  l.E-15 
ELSE 

ADUM  -  A(J) 

BDUM  •  B(J) 

CDUM  -  C(J) 

ENDir 

DENOM  -  1.  -  ADUM*PZII-l) 

PI(I)  -  AOUM/d.  -  AOUM*PZ(  I-l  )  ) 

OZNUM  -  BDUM*P(I,J*l )  ^  CDUM*P( I, J-1  ) 

7  0Z(I)  •  (QZNUM  ♦  ADUM*0I( I-l ) )/DENOM 

DO  8  r  -  MMAX-1,2,-1  (BEGIN  BACKWARDS  SUBSTITION 

DEL  .  PZ( I)*P( I>1,J)  ♦  0Z(I)  -  P(I,J) 
ir(ABS(DCL)  .GT.  ABS ( DELMAX )) THEN 
DELMAX  -  DEL 
IDELMAX  -  I 
JDELMAX  -  J 
ENDir 


64  8  P(1,J)  •  P(I,J)  ♦  RELAX*DEL 

65  10  CONTINUE 

66  RETURN 

67  END 

68 
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E.4  PRESSURE  COEFFICIENT  FORTRAN  PROGRAM 
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Fortran  Program  RESULTS. FOR 

Written  by  Richard  r,  Hubbeil  (July  1987)  Univ.  o£  R.I. 

DIMENSION  R(  101  ) 

COMMON  /r I  ELD/P (  1  250  , ion ,VR(  1250  ,  101  )  ,VZ(  12  50, 101  )  ,  PRESS (  1250.  10/  ) 
data  press. P.VR.VZ, R/1 26 250*0.  ,  1  26250*0  . . 1262 50*0.. 126250*0. .10|*d./ 
DATA  NMA^,.MM>J(,DR,D2,VIN/10l,  12  50.0.00S,  0.02  5, 4  00  ./ 

OPEN(  I, TYPE- 'OLD'  , NAME- ' TDMA . PS  I '  ) 

READ(  1  .  *  )  I  COUNT  j  I  COUNTP  ESMAJC  .  1  MAJ? .  UMAX  ,  RESMAJl 
REAO(  1.8)((P(I,J),  I-l  ,MMAJt)  ,  J-1,  NMAJC) 

8  rORMATI 12625( 10{ ril.8,lX)./) ) 
type  *.'  389.399' 

TYPE  9. ( (P( K,L) ,K.389, 399).L-4,2,-l ) 

TYPE  *. '  400,410' 

TYPE  9.  ( 1 P( K.L)  .E-400,410  ) ,L-4,2.-l ) 

TYPE  *.'  680,690' 

TYPE  9, ( (P( R,L),E-680,690) ,L.4,2,-1) 
type  *,'  691,701' 

TYPE  9, ( ( P{ K,L) ,X.691 . 701 ) /L-4,2,-1 ) 

9  rORnAT(3(lX.llf7.4,/)) 

BOOYMASS  -  0.105/386.1  I LBM-SEC* •2/IN 

BIGKE  -  0.0 

ROW  -  62  ,  4/32.174/20736  .  !WATER  DENSITY  LBM-SEC* • 2/IN* *  4 

TYPE  *,'WHAT  IS  THE  BODY  VELOCITY?  INCHES/SEC' 

ACCEPT  •,VBOD 
DO  3  J  •  l,NnAX 
00  3  I  -  l,MMAX 
3  P(I,J)  -  P( I,J)*VBOD/VIN 

DO  I  J  -  2.NMAX-1 
RtJ)  -  (J-1)*DR 
00  1  I  -  2,MMAJC-l 
tr(i»(I,J)  .GT.  l.C-05)THEN 
VZ(I,J)  -  (P{I.J^i)  -  P{ I , J-l ) )/2 ./R( J )/DR 
VZABS  -  VZ( I, J)  -  VBOD 

vR([,j)  .  -  (p(m.j)  -  p(  :-i ,  J )  )/2 . /R(  J  )/Di 
VELSQ  -  VZ ( I , J ) -VZ (  I , J )  ♦  7R( I , J )  • VR(  I  ,  J  ) 

PRESS(I,J)  -  0.5*ROW*I VB0D**2  -  VELSQ) 

SMALLXE  -  2  .  •  3  . 14 1 59  * R( J ) *0R*02 -ROW* ( VR( 1 , J ) • *2  ♦  V2ABS**2) 

3IGRK  -  BIGKE  *•  SMALLXE 
ENDIf 

1  CONTINUE 

HY?tASS  •  SIGKE/VBOO/VBOD 
TYPE  20,HY?1ASS,BODYMASS 
TYPE  GOING  TO  OUTPUT 

TYPE  *,'  389,399' 

TYPE  9,(  {  P(X.L),K-3a9,399),L-4,2,-n 
TYPE  *,'  400,410' 

TYPE  9,( (P(K,L),K-400,410),L-4,2,-1) 

TYPE  •,'  680,690' 

type  9,( (P(K,L),K.680,690),L-4,2,-1) 

TYPE  *,'  691,701' 

TYPE  9.( (?(K,L)  ,X. 69  1,701 ),L-4, 2,-1 J 

OPEN (2, TYPE- 'NEW' , NAME- ' RESULTS . OUT' ) 

DO  2  I  -  300,300,20 

2  CALL  OUTPUTf  I  , 5,NMAX)  ‘PRINT  OUT  EVERY  flfTH  J  i.t.  J  -  2,7,9 
DO  4  I  -  390,419,1 

4  CALL  OUTPUT! I , 2,NnAX)  ‘PRINT  OUT  EVERY  OTHER  J  !.•.  J  -  2 . 4 , 6  C+C 
DO  5  I  -  420.600.20 

5  CALL  OUTPUT! I , S.NMAX) 

DO  6  I  -  510,698,2 
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6  CALL  OUTPUT! 1,2. NnAX) 

DO  7  I  -  700,800.20 

7  CALL  OUTPUT! I . S.NnAX)  , 

20  rORnAT! ix, 'HYHASS  -  • , f 1 1 . 8 , 3X ,  '  BODYMASS  •  ',fll.8) 
END 


SUBROUTINE  OUTPUT! I . INC . NnAX ) 

COftAON  /r I ELD/P 1 1250, 101), VR (1250, 101), 


V2( 1250,101) , PRESS! 1250, 


HR1TE(2,23)I 


DO  2  JDUH  -  2 ,NnAX. 10*INC 
WRITE!  2,22)  (J,J-0DUn,  9*  INC«-JOUH.  INC) 
HRIT£(2,21)(P( I.J) ,J-JDUH,9*INC^JDUK, INC) 
WRITE!  2,24  )  (VR(  I,J)  ,  J-JDU«,9MNC^JOUH.  INC) 
WRITE(2.24)!VZ( I.J) , J-JDUM.9*INC^JDUH, INC) 
WRITE! 2,21) ( PRESS! I , J ) , J- JOUM , 9 • INC^ JDUH , INC) 
2  CONTINUE 


WRITE! 2, 102) 
RETURN 


21  rORKAT! IX, 10! r7. 4. IX)) 

22  rORnAT(/,lX,10(2X,'J-',l3,lX)) 

23  rORNAT! 30X, ' I  •  ',14) 

24  rORAAT! IX. 10! r7.3. IX) ) 

102  FORMAT!//) 

END 
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